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1. Introduction 

It has been a known fact that laminated fiber composites currently 
in use are relatively weak in resisting impact loads. Great attention 
has been given to modeling the dynamic behavior of composites subjected 
to foreign object impacts and to the search for new forms of composites 
that are capable of improving the impact-resistant property. 

Failure modes in composites resulting from the impacts of a hard 
object and a soft object are, in general, quite different. If the object 
is relatively rigid and small, then the contact time is short and 
extensive damage is usually confined to the neighborhood of the contact 
region.. How to quantify the amount of damage received by the composite 
in the impact zone becomes the central question in the hard-object 
impact problem. 

There are several major factors which could affect the amount of 
damage in a laminated composite due to the impact of a hard object. 

Among them are the mass and approach velocity of the object, the bending 
rigidity of the laminate, and the contact behavior (or the contact law). 
Many researchers have correlated the impact velocity with the damage 
for a given mass. Such relationship between the damage and impact 
velocity becomes invalid if the mass of the striker or the bending 
property of the laminate is changed. The use of a single parameter 
which could account for the combined effect of the above mentioned 
variables is highly desirable. 

Energy dissipation takes place in the process of impact that results 
in damage. It is thus reasonable to use this amount of energy consumed 
in the impact zone to measure the degree of damage in the target compos- 
ite beam. There could be various damage modes such as breakage of 
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fibers, cracking in the matrix, delamination, and plastic deformations, 
which could all contribute to the energy dissipation in the impact zone. 

It is conceivable that analytical estimates of the energies associated 
with these damage modes are prohibitive. A static indentation test 
which produces the loading-unloading curve may prove to be a simple 
means for determining such damage energy, since the energy dissipated 
during the loading and unloading cycle is simply the area enclosed by 
the curve. 

In this report, results of the indentation tests on glass/epoxy 
and graphite/epoxy laminated composites are presented. The results 
show that the loading curve follow the power law with a power index 
1.5, which is identical to the classical Hertzian contact law. Sub- 
stantial permanent deformations are observed even when loaded at very low 
load levels. The unloading curves also follow a power law. 

A high order beam finite element is used for computing the dynamic 
response of laminated composite beams subjected to the impact of an 
elastic sphere. This finite element includes the classical elastic 
Hertzian law of contact as well as the measured contact law. The computer 
program developed for this beam finite element is listed as Appendix A. 

A simple method has been developed for computing the contact force and 
contact duration. An estimate of the contact duration is needed in the 
finite element program in selecting a proper time increment in the time 
integration procedure. This method is found to be quite accurate except 
for very thin beams. 
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2. Indentation Law for Hard Object Impact of Composites 
2.1 Hertzian Law of Contact 

When two solid bodies are in contact, deformation takes place in 
the contact zone and the contact force results. Once the contact force 
is obtained, conventional methods for stress analysis can be used to find 
the stress distribution in the bodies. To determine the contact force - 
indentation relationship often becomes the most important step in analyzing 
the contact problem. 

The most famous elastic contact law was developed by Hertz [1] for the 
contact of two spheres of elastic isotropic materials. The problem was 
solved based on the theory of elasticity. A special case is that if the 
radius of one of the spheres becomes infinite,, then the problem becomes the 
contact of an elastic sphere and an elastic half space. The contact force 
F and the indentation depth a were found to have the relation 


where 




1/2 

s 




( 2 - 1 ) 


( 2 - 2 ) 


In Eq. (2-2), is the radius of the sphere, v is the Poisson's ratio, 

E is the Young's modulus, and the subscripts 1 and 2 indicate the two 
bodies. Equation (2-1) is usually called the Hertzian law of contact for 
a sphere on' half space. 

The 3/2 power law given by Eq. (2-1) was found to be valid by Willis 
[2] for a rigid sphere pressed on a transversely isotropic half space. 

A modified contact law with 
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_ 4 p1/2 
“ 3 




(2-3) 


was employed by Sun [3] for a study on impact of laminated composites. 

In Eq. (2-3), R^, and are the radius, the Poisson's ratio and the 

Young';s modulus of the isotropic shpere, respectively, and Ej is the 
Young's modulus of the fiber- reinforced composite normal to the impact 
plane. 

In applying the classical Hertzian contact law to the impact of 
laminated fibrous composites we face several uncertainties. First, the 
half space assumption is not valid since the laminates in use are of 
finite thickness. Second, the anisotropic and nonhomogeneous property of 
laminated compostes may alter the form of the law. Third, the strain rate 
effect which is not accounted for by the Hertzian law may have significant 
effect on the F-a relation. Except for the strain rate effect, the first 
two uncertainties are soTvable by analyzing the exact contact problem 
of a sphere pressed into a laminated composite using three-dimensional 
elasticity. However, experience tells us that analytical solutions for 
such contact problems are extremely difficult to obtain especially if 
permanent deformations are to be accounted for during unloadings. Since 
unloading paths are particularly important in our study, the experimental 
approach is taken to determine the law of contact for composites. However, 
in this study, the strain rate effect is still neglected. 

2.2 Indentation Law for Laminated Composites 

2.2.1 Theoretical Model 

In this study the general form for the indentation law for laminated 
composites is extended from the classical Hertzian Law. We issume that 
for loading 
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F = k (2-4) 

where k and n will be determined experimentally. It is obvious that when 
n = 3/2 and k is given by Eq. (2-2), this relation becomes the Hertzian 
law for isotropic bodies. It is noted that the constant k has a very 
strange unit if n is not an integer. Also, the value of k depends on 
the unit used for a. A more physically meaningful expression may be 
derived by using a nondimensional indentation depth 

B = a/R^ (2-5) 

with which the indentation can be written as 

F = k* s" (2-6) 


In Eq. (2-6), k* has the unit of force. For the Hertzian law, 



2 

s 




(2-7) 


Permanent indentations in composite targets are usually generated 
even at relatively low projectile impact speeds. The permanent indentation 
accounts for the major part of the energy loss of the projectile. Some 
energy imparted from the projectile to the target can be stored in the 
form of vibrational energy in the target. As far as the local damage 
at the impact zone is concerned, the permanent indentation is of more 
interest to us. For this reason, the force- indentation law for the recovery 
process must be established. In this study, we assume, in the recovery 
process. 




F = F 


m 


«m" “o 


(2-8) 
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where is the maximum contact force just before unloading takes place, 

is the identation corresponding to F^, and is the permanent indentation. 
This recovery law was proposed by Barnhart and Goldsmith [4] for impact 
of a steel ball onto an armor plate. 

2.2.2 Experimental Results 

The experimental set-up is depicted by the sketch in Fig. 2.1. 

The indentation was measured by a dial gage that permits reading up to 
1/5000 in. The dial gage was mounted on the loading piston so that only 
the relative displacement between the indentor and the beam was recorded. 

The indentor was a steel ball of in, diameter. The beam was clamped at 
both ends with various spans. 

Two types of laminated composites have been tested, namely glass/ 
epoxy and graphite/epoxy. The glass/epoxy was Scotch Ply 1002 by the 3M 
Company. It contained 10 0°-plies and 9 90°-plies which alternate in the 
layup with one 0°-ply on top and one at the bottom. The thickness of 
the beam was 0.19 in. and the width was 1,5 in. The graphite/epoxy 
specimens were [0/(j45)2/02/+45]^ laminates. Three different spans, 2 in., 

4 in. and 6 in., were used for the glass/epoxy laminates and two spans, 

2 in. and 4 in., were used for the graphite/epoxy laminates. 

The Loading Curve 

For the glass/epoxy laminate, three sets of loading data were obtained 

for each span. These data were used to determine the best fit for the 
power law, Eq. (2-4), using the least squares method. The results were 

presented in Figs. 2. 2-2. 4. The power indices for the three cases appear to 
be rather close to that of the classical Hertzian law for isotropic media. 
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i.e., n = 1.5. The small deviation from n = 1.5 could be due to measure- 
ment errors. For this reason, we set n = 1.5 and then determined k by 
using the least square fit. The resulting curves are shown in Figs. 

2. 5-2. 7. These curves seem to fit the data very well also. 

The results of the indentation test on the graphite/epoxy laminated 
beams are presented in Figs. 2.8-2.10. For the 2-inch span, the best 
least square fit is n = 1.5; and for the 4-inch span as shown in Fig. 2.10, 
n = 1 .5 also yields a very good fit. 

Table 1 summarizes the indentation laws (the loading portion) obtained 
from the experimental results for a glass/epoxy composite and a graphite/ 
epoxy composite. It is interesting to note that with n = 1.5, the values 
of k for different spans are almost a constant. This indicates that the 
indentation law is ..independent of span. In other words, the bending 
stress does not influence the "contact rig.idity". 

Table 2 presents the indentation laws in terms of B and k* with 
n = 1 .5 (see Eg. (2-6) ). 

The Unloading Curve 

Form the test results we have observed that permanent deformation 
would occur after an indentation test no matter how., small the load was. 

The unloading paths are very different from the loading path as can be 
seen from Figs. 2.11-2.14 for both glass/epoxy and graphite/epoxy. The 
unloading curve is modeled by using Eq. (2-8) in which q and have to 
be determined. Since the permanent indentation depth is difficult to 
measure, the whole data for each unloading path were taken to determine 
the two parameters q and a^. The value q = 2.5 seems to yield the best 
overall fit as shown in Figs. 2.11-2.14. For q = 3.0 (see Figs. 2.15-2.18) 
becomes negative in some cases. 



(a in inches). 


Table 1. Indentation law F = k a’^ 



Glass/Epoxy 

[( 0 / 90 ) 4 / 0 / 90 / 0 /( 90 / 0 )^] 

Graphite/Epoxy 

[0/(145)2/02/145]^ 

Span 

n 

2 

4" 

6" 

2 " 

4" 

Least Squares 
Fit 

n 

1 .54 

1.54 

1 .66 

1.5 

1 .63 

k 

5.569x10^ 

5.603x10^ 

9.655x10^ 

5.964x10® 

9.99x10® 

1 .5 Power 
Fit 

n 

1 .5 

1.5 

1 .5 

1 .5 

1 .5 

k 

4.617x10^ 

4.633x10^ 

4.592x10® 

5.964x10® 

5.126x10® 

Modified 
Hertzian Law 
Eq.(2-3) 

F = 5.461 X 10^ 

Rg = 0.125", Vg = 0.3, Eg = 30 X lO^psi. 
Ej = 1.2 X 10^ psi . 

F =5.24 x 10® a^'® 
Ej = 1.15 x 10®psi. 






















Table 2. Indentation law F = k* 



Glass/Epoxy 

[(0/90)4/0/90/0/(90/0)^] 

Graphite/Epoxy 
[0/( +45) 2 / 02 / 145 ]^ 

Span 

T 1 

2" 

4" 

6" 

2" 

4" 

1 .5 Power 
Fit 

n 

1 .5 

1 .5 

1 .5 

1 .5 

1 .5 

k* 

2.0405x10^ 

2.0475x1 0'^ 

2.0294x10^ 

2.6357x10^ 

2.2654x10^ 

Modified 
Hertzian Law 
Eq. (2-6) 

4 15 

F = 2.4134x10^ g' 

4 15 

F = 2.32 X 10^ 
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Fig. 2.2. 


Least-square fit of the contact force - indentation 
relation for glass/epoxy with 2-inch span. 
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Fig. 2.3. Least-square fit of the contact force-indentation 
relation for glass/epoxy with 4-inch span. 
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Fig. 2.5, 


Least-square fit with n = 1 .5 for glass/epoxy with 
2~inch span. 
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Fig. 2.6. Least-square fit with n = 1.5 for glass/epoxy with 
4- inch span. 
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Fig. 2.7. Least-square fit with n = 1.5 for glass/epoxy with 
6-inch span. 
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Fig. 2.8. Least-square fit of the contact force - indentation 
relation for graphite/epoxy with 2-inch span. 
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Fig. 2.9. Least-square fit of the contact force - indentation 
relation for graphite/epoxy with 4-inch span. 
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Fig. 2.10. Least-square fit with n = 1.5 for graphite/epoxy 
with 4-inch span. 
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Fig. 2.11. Unloading curves for glass/epoxy with 2-inch span. 
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Fig. 2.13. Unloading curves for glass/epoxy with 6-inch span. 
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Fig. 2.17. Unloading curves for glass/epoxy with 6-inch span. 
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3. IMPACT RESPONSES BY FINITE ELEMENT ANALYSIS 
3.1 The Finite Element 

A beam finite element with six degrees of freedom has been developed 
for the dynamic response of elastic isotropic beams subjected to impulsive 
loadings [5], This high order beam element has been shown to be more 
efficient than the conventional element with four degrees of freedom. 

The element displacement function is taken as 

V = a-j + a2X + a^x + a^x + agX + agX (3-1) 


where v is the transverse displacement and a., are constant coefficients. 

The three degrees of freedom at each node are the transverse displacement 
V, the rotation e, and the curvature k. The coefficients a^ in Eq, (3-1) 
can be replaced by the six generalized nodal displacements at the two end 
nodes and, as a result, the displacement function can be alternatively 
expressed in terms of the nodal displacements. 

The stiffness and mass matrices corresponding to the element displacement 
function has been presented elsewhere [5] and are reproduced in the following: 


1200 


[k] 


70L^ 


600L 

30L^ 

-1200 

600L 

-30L^ 

384L^ 

22L^ 

-600L 

21 6l' 



6L^ 

2 

-30L 

8L^ 




1200 

-600L 

P 

30^ 




384L^ 

-22L^ 







(3-2) 
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21720 3732L 28U^ 6000 -181 2L 181 

832L^ 69L^ 1812L -532L^ 52L^ 

4 2 3 4 

6 C 181 r -szr 5 r 

21720 -3732/. 281 

832L^ -69L^ 

6L^ 

where E|^I is the beam bending rigidity, L is the length, p is the mass 
density, and A is the cross-sectional area. If the finite element is 
to be used for the analysis of lamianted composite beams, then the bending 
rigidity Ej^I has to be replaced by the equivalent bending rigidity D. 

3.2 Impact Response 

Based upon the stiffness and mass matrices given by Eqs. (3-2) and (3-3), 
respectively, a computer program has been written specifically for the 
dynamic response of a beam subjected to transverse impact of an elastic 
sphere. A finite difference scheme suggested by Wilson and Clough [6] 
was used to integrate the time variable in the equations of motion. In 
[5], the classical Hertzian law of elastic contact was used to solve a 
few example problems and excellent results were found. 

The finite element program has been modified for the analysis of impact 
of laminated beams. The Hertzian indentation laws, Eqs. (2-1) with Eq. (2-2) 
or Eq. (2-3), as well as the measured indentation formulas can be chosen 
for the analysis. Both elastic loading and actual loading paths can be 
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incorporated in the program. The computer program with a brief user's 
instructions is presented in Appendix A. 

Figures 3.1 - 3.4 show results for some example problems of simply- 
supported steel beams, subjected to impact of a steel ball. The diameter 
of the ball is ^ in. The classical Hertzian law of contact was used in 
the computation. The material constants used are given by Eq. (4-31), 

From these results it can be seen that the impact velocity has a great 
effect on the maximum contact force and contact duration. The thickness 
of the beam has little effect for the two beam depths studied. 

As reported in Section 2, a contact of the steel ball and the glass/ 
epoxy and graphite/epoxy composite always results in a permanent deformation 
The unloading paths are substantially different from the loading path. If 
the actual unloading paths are used, the contact force is certainly expected 
to deviate from that obtained by following elastic unloadings. 

Figures 3.5 and 3,6 present the results for a glass/epoxy laminated 
composite beam with the dimension 0.19 in. D x 1 .0 in. W x 7.5 in. L. 

This is the composite beam used fOr the indentation test. The actual 

5 

indentation law with k = 4.62 x 10 and n = 1 .5 for loading and q = 2.5 
for unloading was used for the computation. Note that, in this case, the 
steel ball has a diameter of 1/4 in. same as that of the identor in the 
static indentation test. The material constants for composite are 

= 5.7 X 10^ psi 
Ej = 1 ,2 X 10^ psi 
Glt " 0-^ ^ 10^ psi 

Vj^T ~ 0.26 

p = 0.002016 s1ug/in^ (0.000168 Ib-sec^/in"^) 
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From the? results in these figures it can be seen that the contact force 
drops more rapidly after reaching its maximum value when the inelastic 
unloading path is followed. However, the total contact duration does 
not seem to be affected by the inelastic unloading. 

The finite element program developed here can also be used in conjunction 
with the experimentally obtained contact law to compute the dynamic strain 
at any point on the beam. The dynamic strain can be experimentally measured 
by using a strain gage. By comparing the measured strain and that predicted 
by the finite element solution, it may be possible for us to determine the 
effect of a result of this comparison. The static indentation law may be 
modified to account for the strain rate effect. The result of the comparison 
will be reported in the future. 
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4. A Simple Method for Computing Contact Force and Duration in 

Elastic Impact 

In using the finite element program described in Section 3, we have 
to choose a proper time increment At and the total length of time integration 
prior to the solution. A poor choice of At may result in poor finite 
difference solutions. A simple way to obtain an approximate impact duration 
prior to the use of the finite element program certainly will avoid futile 
trials. In the following, a simple method is developed for computing 
an approximate contact force and the contact duration. 


4.1 Impact of an Elastic Sphere on a Mass with a Flat Surface 

A simple analysis for a spherical projectile impacting an elastic 
mass with a flat surface was proposed by Timoshenko [7] as follows. 

Denoting the mass and velocity of the target by m^ and v^, respectively, and 
the mass and the velocity of the sphere by m^ and v^, respectively, the 
rates of change of velocity during impact are 


^''t 


(4-1) 


m 



s dt 


F 


(4-2) 


where F is the contact force. The velocity of the relative approach a 
(the indentation) is 



~ V 


t 


From Eqs. (4-1) to (4-3), we obtain 


m. + m 
_t s 

mtm^ 


(4-3) 


a 


F 


(4-4) 
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Substituting the Hertz law of contact, Eq, (2-1), in Eq. (4-4), 
we obtain 


a = - k^a 


3/2 


(4-5) 


where 


K = ~ 


(4-6) 


Integrating Eq. (4-5), we have 


1 ,-2 2^ 2 5/2 

(a - v^) = “ -g- 


(4-7) 


The maximum value of a, a„,„, occurs at a 

MlaX 


0 . We obtain 


max 


( 5 ^ n2/5 
^ 4 kS ^ 


(4-8) 


This together with the Hertzian law yields the maximum contact force. 
From Eq. (4-7), the following relation is derived: 


dt 


da 


(V? - I 


(4-9) 


By introducing 


/4 k c ^2/5 _ a 

(5 2 ' “ - 


(4-10) 


"i 


max 


Equation (4-9) can be rewritten as 


dt = 


max 


dn 


vi (, . „5/2)1/2 


(4-n) 


If we assume that the maximum indentation, is achieved half way 

1TI9X 

through the entire contact, then the duration of impact is obtained from 
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integrating Eq, (4-11) as 


2a 


max 


rl 


'o (1 


dn 


= 2,94 


max 


(4-12) 


4.2 Equivalent Mass Model 

In view of the simple formulas given by Eqs. (4-8) and (4-12), 
we will attempt to find an equivalent mass m^ to represent an actual beam 
or plate. Once this is accomplished, the maximum contact force and the 
contact duration can be estimated easily. 

The equivalent system is developed based upon the condition that 
it stores the same amount of kinetic and strain energies as in the actual 
system. It is assumed that in both systems the strain energies in the 
impactors are negligible and that the kinetic energies are identical. It 
is also assumed that the spheres do the same amount of work on both the 
actual and the equivalent targets. With these assumptions, we conclude 
that the total kinetic energy of the equivalent mass, K^, should be equal 
to the kinetic energy K plus the strain energy U of the actual elastic 
target, i .e. , 


K + U = K, 


(4-13) 


The kinetic energy in the equivalent target system is simply 


h’ht ''I 


(4-14) 


From Eq. (4-1), the velocity of the equivalent mass can be obtained 


by integration as 
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1 

= - — F(t) dr (4-15) 

t •'o 

From all the previous studies, the contact force history resembles 
a sine function. In view of this, we approximate the contact force as 
follows 




(4-16) 


Substituting Eq. (3-16) into Eq. (3-15) and integrating from t=o to t = 
T/2 we obtain the velocity of the equivalent target at t = T/2 as 


V = - i- I F 


(4-17) 


Substitution of Eq. (4-17) into Eq. (4-14) and then into Eq. (4-13) lead to 


(K + U). T/o = i — f"Lv 

t=T/4: 4 m^ IT max 


(4-18) 


Since the deflection of the beam is proportional to the applied force F. 

2 

both U and K contain F terms and can be factored out as 

max 


? ? 

U = F U* K = F K* 
^max ^ ’ ^max ^ 


(4-19) 


in which U* and K* do not depend on F . Equation (3-18) can now be 

max ^ 

written as 

t2 

= (U* + K*). (4-20) 

2m^/ 

From Eqs. (4-6), and (4-8) and (4-12), we note that the contact 
duration T is a function of the equivalent mass m^. Thus, Eq. (4-20) is 
basically a nonlinear equation for m^. Numerical methods will be used to 
find solutions for this equation. 
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4.3 Simply-Supported Beam 

Consider a beam of cross-sectional area A and bending rigidity D. 

The equation of motion is 

“Is + p* ^ <‘'-21) 

3t 

where p is the average mass density (over the thickness) and q(x,t) is a 
time dependent forcing function. For a homogeneous elastic beam, we have 


D = El 


(4-22) 


For laminated composite beams, D is estimated according to Eq. (4-36). 

If the force is a concentrated force F(t) applied at x=c, then the 
solution for Eq. (4-21) can be expressed as [8] 


w(x.|t) 


1 

pA 


CO 


I 

n=l 


Wn(x)Wn(c) 


*■ w^(x)dx 
0 


•t 

F(x)sin w^(t-x)dT 

0 


(4-23) 


where w^(x) is the shape function for the nth natural mode of vibration, 
and is the corresponding natural frequency. 

For a simply-supported beam, we obtain 


w„(x) 


sin 


niTX 

L 


(4-24) 


and 


2 


/niTx4 D 
^ pA 


(4-25) 



43 


If the concentrated force is given by Eq. (4-16), then the beam deflection 
w can be obtained from Eq. (4-23) as 


2 F L. ” 1 

w(x,t) = — I "T 

tt^D n=l n^ 


— (sin Y t 

4t' - 


sinai^t) 


w^(x) , for t£ T 


(4-26) 


In Eq. (4-26), 


T = 


2ir 

“n 


(4-27) 


is the period for the nth mode. The strain energy and the kinetic energy 
can be computed in a straightforward manner. We obtain at t = T/2 


u* . 16 l¥ 

4„ *^1 

IT D 


K* = AlV 

^ 6^2 9] 

7T D 


(4-28) 


where 


n=l 


n 


4n^T^-T^ 


1 - 


sin(n'^ - 5 - ) 


2 

2n T 


w (c) 
n 


(4-29) 


9] = I 

' n =1 


1 


4 2 2 

4n T 


2 ^ 

cos(n — ) w^(c) 


(4-30) 


From the numerical examples, it has been observed that use of fifty 
terms in the series in Eqs. (4-29) and (4-30) should provide a converged 
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solution. In all examples presented in this section, the classical 
Hertzian law is used. 

As a first evaluation of the equivalent mass concept, we consider 
a problem solved by Timoshenko [9] using a numerical procedure to solve 
a nonlinear integral equation. The steel beam considered has a 0.39 in. 

X 0.39 in. (1 cm X 1 cm) cross-section and 6.04 in. (15.35 cm) length. 

The beam is simply-supported at two ends and subjected to impact of a steel 
ball with 0.79 in. (2 cm) diameter. The material properties are 

E = 31 X 10® psi 

V = D.29 (4-31) 

P " 0.00894 s1ug/in^ (0.000745 Ib-sec^/in^) 

It should be pointed out that in the numerical computation, the value for 
the mass density as given by Eq. (4-31) should be divided by a factor of 
12 if the length is given in inches. 

Figure 4.1 shows the contact force histories according to Timoshenko's 
solution and the equivalent mass model. Excellent agreement is noted. 

Figures 4.2 and 4.3 show the contact forces of a simply-supported steel 
beam subjected to impact of a steel ball of 1/2 in. diameter with different 
velocities. The beam has a 1/2 in. x 1/2 in. cross-section and is 30 in. 
long. Both the equivalent mass model results and the finite element results 
are found to have a very close agreement. 

The results for a thicker steel beam (1/2 in. W x 3 in. 0 x 30 in. L) 
with simple supports are presented in Figs. 4.4 and 4.5 for v. = 12 in/sec. and 
1200 in. /sec., respectively. Again, the equivalent mass model works quite 
well in predicting the magnitude and duration of the contact force. 
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Figure 4.6, shows the results for a simple-supported thin steel beam 
(0.5 in. W X 0.08 in. D x 15 in. L) subjected to the impact of a steel ball 
of 0.5 in. diameter with = 100 in. /sec. The equivalent mass model is 
able to predict the maximum contact force but not the contact duration due 
to the long tail portion. 

Figure 4.7 shows the contact force history for a composite beam of the same 
dimension and impact condition as the previous problem. The laminated 
beam consists of 16 piles of graphite/epoxy composite. The ply- thickness 
is 0.005 in. and the lay-up is (0/90/0/90)2^. The material constants are 

= 30 X 10® psi, E.^ = 0.75 X 10® psi 

= 0.4 X 10® psi, = 0.25, (4-32) 

p = 0.00178 slug/in^ (0.000148 Ib-sec^/in^) 

The modified Hertzian law of contact given by Eq. (2-3) was used for the 
solution. Again, from Fig. 4.7 we find that the eqivalent mass model is 
excellent in predicting the maximum contact force but poor in estimating 
the total contact time. From the numerical examples carried out, it seems 
that the equivalent mass model can not yield accurate contact time if the 
target is too thin. 

4.4 Simply-Supported Rectangular Plate 

The plate theory developed by Whitney and Pagano [10] for laminated 
composites is used for the analysis. This plate theory takes the transverse 
shear deformation into account and has been shown by Sun and Lai [11] to 
be adequate for wave propagation. For simplicity, only cross-ply laminated 
plates are considered, for which the equations of motion are given by 
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°n'*'x,xx ^66 '^x,yy ^66^'^y,xy " ^55^x " ^55'^’x " 


(D,,j + DrrJ’P ■•■ Drc 'i' + Roo t - K A«/|\1 j - K A«/,W, = pi't' (4-34) 

' 12 66'^x,xy 66 ^y,xx 22 ^y,yy 44^y 44 ’y ^ ^y ' ^ 


" ''SB'^X.X “ ''55«'xx ^ My.y + ■= ''44’*'yy 

where a comma indicates partial differentiation, q is the lateral load. 


(4-35) 


w is the transverse displacement, i[; and \l) are the rotations of the plane 

A y 

2 

sections, k(=it /12) is a shear correction factor, p is the average mass 

density (over the thickness), h is the plate thickness, I is the rotary inertii 
and 

h/2 


(A. D, .) 


•h/2 


Q^.jd ,z^)dz 


(4-36) 


In Eq, (4-36), are the reduced stiffnesses for the composite material. 
For an isotropic elastic plate, the following relations exist: 

3 


^11 " ^22 


E h' 


12(l-v'') 


Di 2 = vD^., 


D = D 
*^66 2 ^11 


^11 " ^22 


Eh 

1-v' 


(4-37) 


Ai 2 = vA.,,, 


A 


- A - Eh 
44 ~ ^55 " 2(Uv) 


The equations of motion given by Eqs. (4-33) to (4-35) reduce to those for 
the Hindi in's plate theory [12]. 
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If we separate the total displacemert into the bending part, Wj^, and 
that due to the transverse shear deformation, w^, then we have 

♦x = - "b,x 

tij = - w. (4-38) 

b,y 


w = w. + w^ 

D S 


In terms of Wj^ and w,, , the equation of motion can be written as 


^b,xxx ^°66^'^b,xyy ^55^;,x " ^^'^b,x 


(Di 2 + 2Dgg)Wj^^^^y + D22”b,yyy '^44\,y " P^^,y 


^ ^5^,xx ^ ^4^,yy ^ ^ 


Combining equations (4-39) with (4-40), we have 


DllW. „ „ + 2(0,,, + Dc£.)W, + Doo W. + K ArrW 

11 b,xxxx ' 12 66^ b,xxyy 22 b,yyyy 55 s ,xx 


^ ^44'^s,yy *^^^'^b,xx '^b,yy^ 


Equations (4-41) and (4-42) can also be expressed in the form 


2 

L-|Wg + L2Wg ~ pi 5" V w 


3t 


2 ' -"b 


*-2'^s ^ ^ 


9t 


2 "™b ■ "s' 


(4-39) 

(4-40) 

(4-41) 


(4-42) 


(4-43) 



where 


h = “n 3/ * 2(D,2 + 2066) 3^23^2 


+ D 


22 ^ 4 
8y 


^2 ^55 !..2 ^44 ! 2 


3x 


9y 




. 2^2 

ax 3y 


Applying Laplace transform to equations (4-43) yields 


2 2 - 

L-jW, + L,,w = pi s V w. 
lb 2 s b 




(4-44) 


where w^ and w^ are the transformed functions of Wj^ and w^, respectively, 
and s is the Laplace transform parameter. Since the rotatory inertia is 
small, it is neglected in this study. 

Solving Eqs . (4-44), we obtain 


{L 2 " L.| ) ph s'^ + L^L2 


w = (-L, + L„)q 


(4-45) 


We expand the displacement w and the load q in terms of the shape functi 

w„ (x,y) of the natural modes of the plate as 
mn 
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For a simply-supported rectangular plate, the shape function for the (m,n) 
mode is given by 


/ \ . mTrx . niry 


(4-47) 


where a and b are the lateral dimensions of the plate. 
Applying Laplace transform to Eq. (4-46) we obtain 


" = n (s) w^„(x,y) 
m n 


q = n (s) 


(4-48) 


m n 


Substitution of Eqs. (4-48) and (4-47) into Eq. (4-45) leads to 


S + 0) 




(4-49) 


mn 


where 


2 1 ^mn ^mn 


mn h C + E 
mn mn 


= - 2'”i 2 * 2D,,)(e)2(M)2 , 0^^ (-,4 


mtr\2 


'mn 


'55 'a 


(^) 


^ ^4(r)' 


(4-50) 


The quantity is the angular natural frequency for the (m,n) mode. If 
the transverse shear deformation is neglected (i.e. the classical plate theory). 


then 
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U) 


2 _ mn 

mn ~ ph 


(4-51) 


The solution for w can be obtained by applying inverse transform, 
We obtain 


w = ^ II sin 


, , , .in ^2:^^ sin ^ 
ph ^ a b 


rt sin u„„(t-T) 

q (,,) mu dr 

0 mn 


(4-52) 


where 




L. 

ab 


ra 


q(x,y,t) sin sin ^ dxdy 


(4-53) 


Consider a contact force given as a sine function, see Eq. (4-16), 
which is applied at the point (x,| , y,|). Then 

n rb 


\n^^^ ~ ab 


0 


q(x,y,t) sin 2M sip 11^ dxdy 


mirX, 


niry 


4 ^ . ,TTt^ . n 

ib f^max (-t) 


1 


for t < T 


(4-54) 


Substitution of Eq. (4-54) into Eq. (4-52) yields 


4F„ .. CO 00 

w(x,y,t) II 


phab 


m n 


mirx, n-rry, 

(sin sin —r — )(sin sin ^2^) 

a D a D 


1 


■> 

mn 


(.JL_)2 


(sin ^ t - sin w_t) 


mn 


mn 


(4-55) 
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for 0 _< t £ T. For = a/2, = b/2, Eq. (4-55) becomes 


w(x,y,t) = 


4 P 00 00 

max ^ ^ 


phab 


1 


m=l n=l 


(sin sin if) ^ 




T'' 

mn 


2 T ^ ■ Toi) 


mn 


sin oit) 
mn 


mn 


miT niT 

Sin 2 sin ^ 


(4-56) 


comparing Eq. (4-46) with (4-56) we find 


4F , 

R - max 1 _ 

“mn^^^ " phab 2 , 

U) 1 

mn 


j X sin ^ sin x (sin t 

/ TT \2 2 2 T 


mn 


Tto„ 


mn 


sin 


(4-57) 


The kinetic energy in the plate at any time t £ T is given by 


K (t) 


ph 


(•a 


/ 9w \ 2 I I 


(4-58) 


Substituting Eq. (4-57) into Eq. (4-46) and then into Eq. (4-58) we obtain 


K (t) = ^ I I (t) 
8 i: „ mn ' ' 
m n 


(4-59) 


By introducing the stiffness, of the plate system for the (m, n) 
made, the total strain energy can be formally written as 


U (t) = 1 y y K (t) 

' ' 2 t k mn mn ' ' 


m n 


(4-60) 
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Upon substitution of Eqs. (4-59) and (4-60) into the Lagrangian 
equation of motion we obtain 


1 


phab B (t) + K B (t) = Q 

^ mn ' • min mm ' ' 'i 


mn mn mn 


mn 


(4-61) 


where is the generalized force. From Eq. (4-61) we obtain the natural 
frequency for the (m,n) made as 


K 


mn phab mn 


(4-62) 


from which 


mn 


phab 2 
4 “mn 


(4-63) 


Substituting Eqs. (4-63) and (4-57) into Eq. (4-60) we obtain 




1 1 

m n 


1 


1 


"mn 1 - (-^) 


. miT . niT 

o X Sin — o sin -o X 

TT \ Z 2 L 


X (sin ^ t - V- sin w t) 

T mn 


(4-64) 


With Eqs. (4-60) and (4-64), the quantities U* and K* in the equivalent 


mass mode can be obtained. We have 
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U* (T/2) = 


K* (T/2) = 


phab 3 


ilL 
phabT 


2 ^3 


(4-65) 

(4-66) 


where 


f, = n 


m n 


1 1 


"mn 1 - 


-JL— V 
mn 


(1 - 


Tw, 


sin -?^) sin 


mn 


mir mr 
2 s T n 2 


(4-67) 


93 = n 
m n 


1 1 


1 - (^) 
mn 


ui T 

cos (-f^-) sin ^ sin ^ 


(4-68) 


Karas [13] considered the impact of a steel ball of 2 cm in diameter 
on a simply-supported square steel plate with a=b=20 cm and h = 0.8 cm by 
using the classical plate theory. The impact velocity of the ball was 100 
cm/sec. The contact force histories obtained by Karas and by using the 
equivalent mass model are shown in Fig. 4.8. It is seen that the equivalent 
mass model yields a good estimate of the maximum contact force and contact 
duration. 
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Fig. 4,6 Simply-supported steel beam (0.5"W x 0.08"D 
X 15"L) subjected to impact of a steel ball 
at 100 in/sec. 
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5. Conclusions 

Static indentation tests have been performed to determine the law 
of contact between a steel ball and two laminated composites, namely, 
glass/epoxy and graphite/epoxy. It has been found that the loading 
path followed very well the power law 

p - ^ 1-5 

r - k a 

where F is the contact force, k is a coefficient, and a is the indentation 
depth. Tests were conducted with beams clamped at two ends with various 
spans. The results indicated that the indentation law does not seem to 
depend on the span between the clamps. The experimental results have 
also revealed that both composites tested possessed a pronounced inelastic 
behavior even at very low contact force levels. The unloading paths 
from various loading points have been obtained experimentally and fitted 
into a power law for the computational purpose. 

An efficient high order beam finite element has been employed together 
with the classical Hertzian contact law or the measured contact law for 
analyzing the impact response. The finite element program is capable 
of computing the contact force, contact duration, and all the dynamic 
responses in the laminated composite. A simple method for estimating the 
contact force and duration has been developed and shown to be quite 
accurate except for very thin beams. 
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APPENDIX A 

A COMPUTER PROGRAM FOR FINITE ELEMENT ANALYSIS OF THE TRANSVERSE IMPACT 
OF A BEAM 

The following is a description of the input data required to 
analyze the transverse impact of a beam. The description is by card 
sections, and where applicable, the number of cards precedes the name. 
The arrangement of the cards is shown in Fig. A-1. 

Heading Card (s) (12, 10A7) 

One card is required for each problem. 

Cols. 1-2 Problem number (NPROB) 

3-72 Arbitrary problem identification (TITLE) 

2. 1-Control Card (915) 

Cols. 1-5 Number of nodal points (NP) 

6-10 Number of elements (NE) 

11-15 Number of restrained boundary nodes (NB) 

16-20 Number of output printing cycles (NTM) 

21-25 Number of material types (NMAT) 

For isotropic materials, this number is limited to 24 
plus one for the sphere. However, for a laminated 
composite, this number can only be two. 

26-30 Output printing frequency in ysec (NDIN) 

31-35 Beam material type (MATP) 

0 - if beam is isotropic 

1 - if beam is a laminated composite 
36-50 Number of nodal data cards (NDC) 

Explained later. 
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41-45 Control for print of input data (11) 

0 - Input printed at beginning of first problem only. 

1 - Input printed for each new problem. 

Input for the printing scheme outlines the cycle and frequency at 
which the output is printed. The integer, NTM, indicates how many 
times output is printed after the sphere makes contact and the integer, 
NDIN, indicates how much time elapses between printing of the output. 

In addition, NDIM is measured in tenths of a microsecond. As an example, 
if one wishes to print output every 5 ysec for 10 cycles, then NDIN 
equals 50 (in ysec) and NTM=10. Observe that (NDON x NTRM)/10 yields 
the time at which computations stop, in this case its 50 ysec. 

3 . 1 - Dimension Card ( 3 FI 0 . 0 ) 

Cols. 1-10 Beam thickness (TB) 

1 1-20 Beam width (WB) 

21-30 Sphere radius (R) 

4. 1 - Nodal Impact Card (I5,2F10.0) 

Cols. 1-5 Impacted node (NQ) 

6-15 Impact velocity (Q2) 

16-25 Time increment (DT) 

5. Element Type Material Properties Card (s) (I5,5F10.0) 

1 card per material 

Cols. 1-5 Material number (IMAT) 

6-15 Longitudinal Young's modulus (0RT(N,1)) 

16-25 Transverse Young's modulus (0RT(N,2)) 

26-35 Shear modulus (0RT(N,3)) 

36-45 Poisson's ratio (0RT(N,4)) 

46-55 Mass density, p (0RT(N,5)) 
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The last material card must contain the material properties of 
the impacting sphere. If the sphere and the beam possess identical 
material properties, then only one material card (NMAT = 1) is necessary. 

6. 1 - Identation Law Card (E10.3,2F10.0) 

Cols. 1-10 Loading coefficient k (STF) 

11-20 Permanent deformation (DISPEM) 

21-30 Unloading power q (QP) 

If the Hertzian law is used for loading, set STF = 0.0. If elastic 
unloading is followed, then set DISPEM = 0.0 and the input for QP will 
be ignored. 

7. Nodal Data Card(s) (2l5, 2F10.0, I5) 

1 card is required for each set of identical elements. 

Cols. 1-5 Beginning node in the set (NDl) 

6-10 Final node in the set (ND2) 

11-20 x-position of beginning node (XI) 

21-30 x-position of final node (X2) 

31-35 Element material type of set (IMT) 

This input provides information for the automatic element generator 
in the program. Given the above information for each set of identical 
elements, the program computes the x-position of each node and assigns 
each element a material type and the Ith and Jth nodes. The number of 
these cards is equal to NDC, which is input on the control card. 

NOTE: Node 1 must begin at position x - 0. 

8. Boundary Conditions Card(s ) (215) 

1 card per restrained node 

Cols. 1-5 Restrained node (NBC) 

6-10 Boundary condition code (NFIX) 

The boundary condition code is an integer containi r.g three digits. 



67 


Each digit in the code is either 1-restrained or 0-free, The ones 
digit controls the curvature, the tens digit controls the slope, and 
the hundreds digit controls the displacement. As an example, if one 
node was clamped, then the displacement and slope are zero and the 
curvature is nonzero, or 

v = 0 0 = 0 therefore Code 110 

NOTE; Boundary conditions may be specified at any node with any code. 

8. Number of layers in laminate (15) (MLAYER) 

9. Laminate data (15, F5.0, FlO.O) 

1 card per layer. 

Cols. 1-5 Layer number (L) 

6-10 Fiber orientation (TH) 

11-20 Layer thickness (TK) 

If a laminated composite beam is to be examined under impact, two 
major alterations in the program must be made. This program provides 
for both, with the proper indication on the control card (MATP =1). 

From thg laminate data given, an equivalent bending rigidity is computed, 
or D-|i = El. In addition, the contact coefficient in the Hertzian 
Contact Law is also computed differently for composite beams. 

NOTE: If an isotropic beam is used, skip Cards 8 and 9. 

10. Termination Card 
EXAMPLE 1 

Consider the impact of a steel sphere on a steel cantilever beam. 

The dimensions of the beam are 0.5" W x 0.08" D x 15" L and the sphere 
has a diameter of 0.5" in. The initial velocity of the sphere is 100 
in/sec., with the point of impact located at the mid-point. Numerical 
solutions are to be obtained up to 100 ysec by using 30 finite elements. 
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The material constants used in this computation are 

E = 30 X lO^psi, V = 0.25 and p = 0.00880 slug/in^(0. 000733 lb - sec^/in^). 
Note that the value of p in slug should be divided by 12 if length is 
measured in inches. 

The sample inputs and outputs for Example 1 and Example 2 are 
listed following Fig. A-1, The results for the contact force, the dis- 
placement of the sphere and the deflection of the beam at the impact point 
are shown in Fig. A-2. The displacement profiles of the beam at the 
impact point are shown in Fig. A-2. The displacement profiles of the 
beam at various times are presented in Fig, A-3. 

EXAMPLE 2 

This example is identical to the previous example except that the 
beam is now a laminated composite which consists of 16 layers of 
graphite/epoxy composite. The ply-thickness is 0.005" and the lay-up 
is ( 0 / 90 / 0 / 90 ) 2 ^. The material constants are 

E.|.| = 30 X 10^ psi E 22 = 0.75 x 10^ psi 6.,2 = 0.4 x 10^ psi 

v ^2 = 0-25 , p = 0.00178 slug/in^(0, 000148 lb - sec^/in^) 

The corresponding results are shown in Figs. A-4 and A-5. 
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Fig, A-1 Deck set-up 





70 


Sample Input for Example 1 











Sample Output for Example 1 


1 0.08DX0.5UX15L ISO. CAMTILEUER WITH 02=100 IM/SEC. 


MODAL f='OIMTS 31 

ELEMEMTS 30 

BOUMDARY CONDITIOMS 1 

OUTPUT LIMIT 1000 

DEGREES OF FREEDOM 3 

MATERIALS 1 


BEAM THICKNESS .080 

BEAM WIDTH .500 

SPHERE DENSITY .000733 

SPHERE RADIUS .250 


IMPACT NODE IG 

IMPACT UELOCITY 100.0 

INTEGRATION TIME INCREMENT! X E-OS SEC) 3.500E-08 
MATERIAL PROPERTIES 


MAT. NO. 

El 

E2 

G12 

U12 

RHO 

1 

30000000.0 

30000000.0 

11500000.0 

.300 

.000733 


PERMANENT DEFORMATION(IN) *.000000 
NODAL POINTS 



X 

V 

1 

0.000 

0.000 

2 

.500 

0.000 

3 

1.000 

0.000 

4 

1.500 

0.000 

5 

2.000 

0.000 

6 

2.500 

0.000 

7 

3.000 

0.000 

8 

3.500 

0.000 

3 

4.000 

0.000 

10 

4.500 

0.000 

11 

5.000 

0.000 

12 

5.500 

0.000 

13 

G.OOO 

0.000 

14 

G.500 

0.000 

15 

7.000 

0.000 

IG 

7.500 

0.000 

17 

8.000 

0.000 

18 

8.500 

0.000 

IS 

9.000 

0.000 

20 

3.500 

0.000 

21 

10.000 

0.000 

22 

10.500 

0.000 

23 

11.000 

0.000 

24 

11.500 

0.000 

25 

12.000 

0.000 

EG 

12.500 

0.000 

27 

13.000 

0.000 

23 

13.500 

0.000 

29 

14.000 

0.000 

30 

14,500 

0.000 

31 

15.000 

0.000 



ELEMENTS 



I 

J 

K 


MAT 

1 

1 

2 

0 

0 

1 

2 

2 

3 

0 

0 

1 

3 

3 

4 

0 

0 

1 

4 

4 

5 

0 

0 

1 

5 

5 

6 

0 

0 

1 

6 

6 

7 

0 

0 

1 

7 

7 

8 

0 

0 

1 

8 

8 

9 

0 

0 

1 

9 

9 

10 

0 

0 

1 

10 

10 

11 

0 

0 

1 

11 

11 

12 

0 

0 

1 

12 

12 

13 

0 

0 

1 

13 

13 

14 

0 

0 

1 

14 

14 

15 

0 

0 

1 

15 

15 

IG 

0 

0 

1 

IG 

IG 

17 

0 

0 

1 

17 

17 

18 

0 

0 

1 

18 

18 

19 

0 

0 

1 

19 

19 

20 

0 

0 

1 

20 

20 

21 

0 

0 

1 

21 

21 

22 

0 

0 

1 

22 

22 

23 

0 

0 

1 

23 

23 

24 

0 

0 

1 

24 

24 

25 

0 

0 

1 

25 

25 

2G 

0 

0 

1 

26 

2S 

27 

0 

0 

1 

27 

27 

28 

0 

0 

1 

28 

28 

29 

0 

0 

1 

29 

29 

30 

0 

0 

1 

30 

30 

31 

0 

0 

1 


BOUNDARY CONDITIONS 

31 no 


PRINTING SCHEME 

1. REPORT OUTPUT EUERY 5.00 MSEC 

2. TERMINATE OUTPUT AT ID 0.00 MSEC 


TYPICAL STIFNESS MATRIX OF AN ELEMENT 


8.777E+04 
2. 134E+04 
5.4SSE+02 
-8.777E+04 
2. 134E+04 
-5.48GE+02 


2, 194E+04 
7.022E+03 
2,011E+02 
-2.194E+04 
3.950E+03 
-7.314E+01 


5.48GE+02 

2.011E+02 

2.743E+01 

-5.48SE+02 

7.314E+01 

4.571E+00 


-8.777E+04 
-2. 194E+04 
-5.48GE+02 
8.777E+04 
-2. 194E+04 
5.48SE+02 


2. 194E+04 
3.950E+03 
7.314E+01 
-2. 194E+04 
7.022E+03 
-2,011E+02 


-5.48GE+02 

-7.314E+01 

4.571E+00 

5.4S6E+02 

-2.011E+02 

2.743E+01 


TYPICAL MASS MATRIX OF AN ELEMENT 


5.743E~08 
4.934E-07 
1 . 85SE~03 
1.587E--0G 


4.934E-07 

5.500E-03 

2.281E-09 

2.39SE-07 


l,858E-0a 

2.281E-09 

9.31GE-11 

l,197E-08 


1.587E-0S -2.39GE-07 
2.38SE-07 -3.517E-03 
1.197E-08 -1.71SE-03 
5.743E-0S -4.934E-07 


-2.3SSE-07 

1.1S7E-03 


-3.517E-08 

l,719E-09 


-1.719E-09 

8.2G3E-11 


-4.934E-07 

1.858E-08 


5.500E-08 

-2.281E-09 


l.lS7E-08 

1.71SE-09 

8.2S3E-11 

1.85SE-08 

-2.281E-09 

9.91GE--11 



73 


0.08DX0.5WX15L ISO. CANTILEUER WITH 02=100 IN/SEC. 


TIME 

ELAPSED(MSEC) 

10.500 



FORCE (LB) 


1.684E+02 



MASS DISPLACEMENT(IN) 

9.747E-04 



MASS 

UELOCITY(IN/SEC) 

7.849E+01 



MASS 

ACCEL 

.(IN/SEC2) 

-3.509E+0B 



INDEMTATION(IN) 

6.197E-04 



NODE 


DISP 

STRAIN-XX 

STRAIN-VY 

STRESS-XX 

1 


1.670E--11 

-3.892E-11 

2.3B8E-11 . 

-2.3B8E-03 

2 


-2.051E--11 

5.506E-10 

-1.B52E-10 

1.B52E-02 

3 


3.351E-11 

-9.258E-10 

2.777E-10 

-2.777E-02 

4 


-3.684E-11 

3.027E-10 

-9.081E-11 

9.081E-03 

5 


-1.222E--10 

5.425E-03 

-1.B28E-0S 

1.628E-01 

6 


8.400E-10 

-2.698E-08 

8.093E-0S 

-8.0S3E-01 

7 


-2.888E-09 

8.059E-08 

-2.418E-08 

2.418E+00 

8 


6.521E-0S 

-1.550E-07 

4.B50E-08 

-4.B50E+00 

S 


-6.831E--0S 

8.597E-08 

-2.579E-08 

2.573E+00 

10 


-1.677E--08 

6.471E-07 

-1.941E-07 

1.941E+01 

11 


1.002E--07 

-2.525E-06 

7.574E-07 

-7.574E+01 

12 


-1.803E--07 

2.561E-06 

-7.B83E-07 

7.BB3E+01 

13 


-4.412E-07 

1.034E-05 

-3.101E-0B 

3.101E+02 

14 


-3.972E-07 

2.215E-05 

-B.B45E-0B 

B.B45E+02 

15 


-3.358E--05 

1.188E-04 

-3.5B5E-05 

3.5B5E+03 

16 


3.578E-04 

-7.110E-04 

2.133E-04 

-2.133E+04 

17 


-3.358E--05 

1.188E-04 

-3.565E-05 

3.5B5E+03 

18 


-3.972E-07 

2.215E-05 

-B.B45E-0B 

B.B45E+02 

IS 


-4.412E--07 

1.034E-05 

-3.101E-06 

3.101E+02 

20 


-1.803E-07 

2.5B1E-0B 

-7.B83E-07 

7.BB3E+01 

21 


1.002E--07 

“2.S25E-06 

7.574E-07 

-7.574E+01 

22 


-1.677E--08 

B.471E-07 

-1.941E-07 

1.941E+01 

23 


-6.891E-0S 

8.597E-08 

-e.579E-08 

2.57SE+00 

24 


6.521E-0S 

-1.550E-07 

4.B50E-08 

-4.B50E+00 

25 


-2.888E--0S 

8.059E-08 

-2.418E-08 

2.418E+00 

26 


8.400E--10 

-2.698E-08 

8.093E-OS 

-8.093E-01 

27 


-1.221E-10 

5.424E-0S 

-1.B27E-09 

1.B27E-01 

28 


-3.6S7E--11 

3.067E-10 

-3.201E-11 

9.201E-03 

2S 


3.985E--11 

-9.352E-10 

2.80BE-10 

-E.80BE-02 

30 


-2.116E-11 

5.B98E-10 

-1.709E-10 

1.709E-02 

31 


1.306E--21 

-2.755E-10 

8.HB4E-11 

-8.2B4E-03 
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0.08DXO.SMX15L ISO. CANTILEUER WITH 02=100 IN/SEC. 


TIME 

ELAPSED (MSEC) 

35-000 



FORCE (LB) 


3.223E-01 



MASS 

DISPLACEMENT(IN) 

2.332E-03 



MASS 

UELOCITYC IN/SEC) 

4.839E+01 



MASS 

ACCEL 

.(IN/SEC2) 

-G.718E+03 



INDENTATION(IN) 

1 . 139E-05 



NODE 


DISP 

STRAIN-XX 

STRAIN-YY 

STRESS-XX 

1 


~l.B39E-08 

1.939E-07 

-5.817E-08 

5.817E+00 

2 


“4.074E-08 

1.315E-0B 

-3.944E-07 

3.944E+01 

3 


1.731E-07 

-3.1B3E-0B 

9.490E-07 

-3.490E+01 

4 


8. 057E-08 

-1.543E-0B 

4,B30E-07 

-4.B30E+01 

5 


-1.G23E-07 

1.457E-0B 

-4.371E-07 

4.371E+01 

6 


“5.1G1E-07 

G.35BE-0B 

-1.907E-0B 

1.907E+02 

7 


-8.232E-07 

8.349E-0B 

"2.505E-0B 

2.505E+02 

8 


-1.B21E-06 

1.321E-05 

-3.964E-0B 

3.9B4E+02 

9 


-3.188E-0G 

1.493E-05 

-4.478E-0B 

4.478E+02 

10 


~l.B13E-07 

-1.B02E-05 

4.806E-0G 

-4.80BE+02 

11 


2.248E-05 

-5.990E-05 

1.797E-05 

-1.797E+03 

12 


-4.545E-05 

1.443E-04 

-4.329E-05 

4.329E+03 

13 


3.28GE-05 

-2.545E-04 

7.63GE-05 

-7.G3BE+03 

14 


“3.94GE-04 

4.B3BE-04 

~1.391E-04 

1.391E+04 

15 


1.079E-03 

“B.524E-0S 

1.957E-05 

-1.957E+03 

IG 


2.322E-03 

-4.584E-04 

1.375E-04 

-1.375E+04 

17 


1.079E-03 

-G.524E-05 

1.957E-05 

-1.957E+03 

18 


~3.94BE~04 

4.B3BE-04 

~l-391E-04 

1.391E+04 

19 


9.28BE-05 

-2.545E-04 

7,S3BE-05 

-7,B3BE+03 

20 


-4.545E-05 

1.443E-04 

"4. 3295-05 

4.329E+03 

21 


2.248E-05 

-5.990E“05 

1.797E-05 

-1.797E+03 

22 


-1.B12E--07 

-1.G02E-05 

4.807E-0G 

-4.807E+02 

23 


-3. 188E-0B 

1.493E-05 

-4.478E-0G 

4.478E+02 

24 


-1.B21E-0B 

1.322E-05 

“3.9GGE-0B 

3.9BBE+02 

25 


-8.225E-07 

8.333E-0B 

-2.500E-0G 

2.500E+02 

26 


-5.171E-07 

G.379E~0B 

-1.914E-0G 

1.914E+02 

27 


-1.B20E-07 

1.454E-0G 

“4.3G2E-07 

4.3B2E+01 

28 


8.2B3E-08 

-1.598E-0B 

4.793E-07 

-4.733E+01 

29 


1.B84E-07 

-3.0S7E-0B 

9.171E-07 

“9. 171E+01 

30 


-3.759E-08 

1.254E-0B 

-3.7G2E-07 

3.7B2E+01 

31 


-2. 132E-18 

5.192E-08 

-1.558E-08 

1.558E+00 
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Sample Input for Example 2 



































76 


Sample Output for Example 2 


1 0.08DX0.5WX15L COMP, CANTILEUER WITH 02=100 IN/SEC 


NODAL POINTS 31 

ELEMENTS 30 

BOUNDARY CONDITIONS 1 

OUTPUT LIMIT 1000 

DEGREES OF FREEDOM 3 

MATERIALS 2 


BEAM THICKNESS .080 

BEAM WIDTH- .500 

SPHERE DENSITY .000733 

SPHERE RADIUS; .250 


IMPACT NODE 16 

IMPACT UELOCITY 100,0 

INTEGRATION TIME INCREMENT! X E-OB SEC) l.OOOE-07 


MATERIAL PROPERTIES 


MAT. NO. 

El 

E2 

G12 

UlS 

RHO 

1 

30000000.0 

750000.0 

400000.0 

.250 

.000148 


2 30000000.0 30000000.0 11500000.0 ,300 .000733 


PERMANENT DEFORMATIONdN) *.000000 


ABD MATRIX 


1.832E+0B 

1.502E+04 

~5.912E-33 

1.502E+04 

l,232E-<-0G 

-1.994E-09 

-5.912E-39 

-1.394E-09 

3.200E+04 

-2.910E-11 

3.183E-12 

-3.788E-55 

3.183E-12 

3,402E-10 

-4.B53E-25 

-9.788E-55 

-4.G53E-25 

4.547E-12 

NODAL POINTS 

X 

1 0.000 

V 

0.000 

2 

,500 

0.000 

3 

1.000 

0.000 

4 

1.500 

0.000 

5 

2.000 

0.000 

6 

2.500 

0.000 

7 

3.000 

0.000 

8 

3.500 

0.000 

9 

4.000 

0.000 

10 

4.500 

0.000 

11 

5.000 

0.000 

12 

5.S00 

0.000 

13 

B.OOO 

0.000 

14 

6.500 

0.000 

15 

7.000 

0.000 

16 

7.500 

0.000 

17 . 

8.000 

0.000 

18 

8.500 

0.000 

19 

3.000 

0.000 

20 

9.500 

0.000 

21 

10.000 

0,000 

22 

10.500 

0.000 

23 

11.000 

0.000 

24 

11.500 

0.000 

25 

12.000 

0,000 

2B 

12.500 

0.000 

27 

13.000 

0.000 

28 

13.500 

0,000 

29 

14.000 

0.000 

30 

14.500 

0.000 

31 

15.000 

0.000 


-2.910E-11 3.183E-12 -9.788E-55 

3.183E-12 3.402E-10 -4.653E-25 

-9.788E-55 -4.653E-25 4.547E-1H 

7.742E+02 8.013E+00 -2.5G2E-42 

8.013E+00 5.338E+02 -8.B38E-13 

-2.5G2E-42 -8.S39E-13 1.707E+01 



ELEMENTS 



I 

J 

K 


MAT 

1 

1 

2 

0 

0 

1 

2 

2 

3 

0 

0 

1 

3 

3 

4 

0 

0 

1 

4 

4 

5 

0 

0 

1 

5 

5 

6 

0 

0 

1 

6 

B 

7 

0 

0 

1 

7 

7 

8 

0 

0 

1 

8 

8 

9 

0 

0 

1 

9 

9 

10 

0 

0 

1 

10 

10 

11 

0 

0 

1 

11 

11 

12 

0 

0 

1 

12 

12 

13 

0 

0 

1 

13 

13 

14 

0 

0 

1 

14 

14 

15 

0 

0 

1 

15' 

15 

16 

0 

0 

1 

16 

16 

17 

0 

0 

1 

17 

17 

18 

0 

0 

1 

18 

18 

19 

0 

0 

1 

19 

19 

20 

0 

0 

1 

20 

20 

21 

0 

0 

1 

21 

21 

22 

0 

0 

1 

22 

22 

23 

0 

0 

1 

23 

23 

24 

0 

0 

1 

24 

24 

25 

0 

0 

1 

25 

25 

26 

0 

0 

1 

2B 

26 

27 

0 

0 

1 

27 

27 

28 

0 

0 

1 

28 

28 

29 

0 

0 

1 

29 

29 

30 

0 

0 

1 

30 30 

BOUNDARY 

31 110 

31 0 

CONDITIONS 

0 

1 


PRINTING SCHEME 

1. REPORT OUTPUT EUERY 5.00 MSEC 

2. TERMINftTE OUTPUT ATIOO.OO MSEC 

TYPICAL STIFNESS MATRIX OF AN ELEMENT 
1.0B2E+05 2.B54E+04 B.B3BE+02 -1.0B2E+05 2.B54E+04 -6 

2.B54E+04 8.434E+03 2.433E+02 -2.B54E+04 4.F78E+03 -8 

G.B3BE+02 2.433E+02 3.318E+01 -B.B3BE+02 8.848E+01 5 

-1.0B2E+05 -2.B54E+04 -B.B3BE+02 1.0B2E+05 -2.B54E+04 S 

a.B54E+04 4.778E+03 8.848E+01 -2.654E+04 8.4S4E+03 -2 

-B.63BE+02 -8.848E+01 5.530E+00 6.B3BE+02 -2.433E+02 3 

TYPICAL MASS MATRIX OF AN ELEMENT 
I.IBOE-OB 9.9B3E-08 3.751E-09 3.203E-07 -4.837E-08 2 

9.963E-08 l.lllE-08 4.B05E-10 4.837E-08 -7.101E-09 3 

3.751E-09 4.B05E-1C) 2.002E-11 2.418E-09 -3.470E-10 1 

3.203E-07 ■ 4.837E-08 2.41GE-03 I.IGOE-OB -9.9G3E-08 3 

-4.837E-08 -7.101E-09 -3.470E-10 -9.963E-08 1.111E-08 -4 

2.41GE-09 3.470E-10 1.BS8E-11 3.751E-0S -4.B05E-10 2 


B3BE+D2 

848E+01 

530E+00 

B3BE+02 

433E+02 

318E+01 


416E-09 

470E-10 

BB8E-11 

751E-09 

B05E-10 

002E-11 



0.08DX0.5WX15L COMP. CANTILEUER WITH 02=100 IM/SEC 


TIME 

ELAPSED (MSEC) 

20.000 



FORCE (LB) 


2.989E+01 



MASS 

DISPLACEMENT(IN) 

1.9B3E-03 



MASS 

UELOCITY( IN/SEC) 

9.407E+01 



MASS 

ACCEL 

.( IN/'SEC2) _ 

-6,230E+05 



INDENTATION(IN) 

1.S57E-03 



NODE 


DISP 

STRAIN-XX 

BTRAIN-YY 

5TRESS-XX 

1 


4.354E-03 

-1.203E-08 

3.008E-09 

-•3.B10E-01 

2 


1.211E-09 

-1.140E-08 

2.850E-09 

-3.420E-01 

3 


1.5S4E-03 

-1.941E-08 

4.853E-03 

-5.824E-01 

4 


4.0B3E-03 

-5.353E-08 

1.338E-08 

-1.60BE+00 

5 


1.002E-08 

-1.048E“07 

2.B19E-08' 

-3.143E+00 

6 


8.16SE-08 

-1.112E-07 

2.780E-08 

“3.336E+00 

7 


2.325E-08 

“4.028E-08 

1.007E-08 

-1.208E+00 

8 


-7.712E-08 

3.099E-07 

-7.748E-08 

9.237E+00 

9 


-2.B02E-08 

-1.B33E-07 

4.081E-08 

-4.898E+00 

10 


3,BB6E-07 

"2.982E-07 

7.454E-08 

-8.345E+00 

11 


~1.131E“0B 

8.802E-07 

-2.201E-07 

2.B4IE+01 

12 


2.941E-0S 

-8.390E"08 

2.098E-08 

-2.517E+00 

13 


-1.8B2E-0B 

~3.902E”Oe 

2.476E-0B 

“2.971E+02 

14 


-3.526E-05 

3.434E-05 

~8.58BE-0B 

1.030E+03 

15 


1.289E-04 

3.957E-05 

-9.892E-06 

1.187E+03 

IB 


4.115E-04 

-1.934E-04 

4.835E-05 

-5.802E+03 

17 


1.283E-04 

3.9S7E-05 

“3.892E-0B 

1.187E+03 

18 


-3.52BE-05 

3.434E-05 

-8.586E-0B 

1.030E+03 

19 


“l.a6e£"0S 

-9.302E~0B 

2.47BE-06 

-2.971E+02 

20 


2.941E-06 

“8.370E--08 

2.092E-08 

-2.511E+00 

21 


-1.131E-0B 

8.798E-07 

-2.199E~07 

2.B33E+01 

22 


3.B6SE--07 

“2,978E-07 

7.44BE-08 

-8.335E+00 

23 


-2.S03E-08 

-l.S23E“07 

4.073E-08 

-4.888E+00 

24 


-7.707E-08 

3.08BE-07 

“7.721E-08 

9.2B5E+00 

25 


2.321E-08 

“3.940E-08 

9.850E-09 

-1.182E+00 

26 


H.1B4E~08 

“l,112E-07 

2.779E-08 

-3.335E+00 

27 


1.003E-08 

“1.050E-07 

2.624E-08 

-3.149E+00 

28 • 


4.071E-09 

-5.328E-08 

1.332E-08 

“1.599E+00 

29 


1.B96E-09 

-2.1B2E-08 

5.40BE-09 

-B.487E“01 

30 


1.136E-09 

-9.681E-09 

2.420E-09 

-2.304E-01 

31 


2.507E-18 

“3.823E-08 

9.5SBE-09 

-1.147E+00 
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0.08DX0.5WX15L COMP. CANTILEUER WITH 03=100 IN/SEC 


TIME ELAPSED(MSEC) 
FORCE (LB) 

MASS DISPLACEMENT(IN) 
MASS UELOCITYdN/SEC) 
MASS ACCEL.(IN/SEC2) 
IMDEMTATIOMCIN) 


30.000 

2.785E+00 

7.288E-03 

6.770E+01 

-5.805E+04 

3.264E-04 


NODE 

DISP 

STRAIN-XX 

STRAIN-YY 

STRESS-XX 

1 

3.708E-05 

-2.198E-08 

5.496E-09 

-6.595E-01 

2 

-3.087E-05 

2.453E-05 

-6.132E-06 

7.358E+02 

3 

2.007E-05 

-1.337E-06 

4.844E-07, 

-5.812E+01 

4 

5.413E-05 

-3.431E-05 

8.577E-06 

-1.029E+03 

5 

-7.656E-05 

1.602E-05 

-4.004E-06 

4.805E+02 

6 

-1.137E-04 

5.467E-05 

-1.367E-05 

1.640E+03 

7 

1.365E-04 

-3.582E-06 

8.954E-07 

-1.074E+02 

8 

3.434E-04 

-9.696E-05 

2.424E-05 

-2.909E+03 

9 

5.487E-06 

-7.98SE-05 

1.996E-05 

-2.396E+03 

10 

-7.572E-04 

6.927E-05 

-1.732E-05 

2.078E+03 

11 

-l.lOOE-03 

1.993E-04 

-4.982E-05 

5.978E+03 

12 

-2.664E-04 

1.998E-04 

-4.995E-05 

5.994E+03 

13 

1.748E-03 

7.519E-05 

-1.880E-05 

2.256E+03 

14 

4.217E-03 

-8.072E-05 

2.018E-05 

-2.422E+03 

15 

6.202E-03 

-1.939E-04 

4.998E-05 

-5.998E+03 

16 

6.969E-03 

-2.614E-04 

6.536E-05 

-7.843E+03 

17 

6.202E-03 

-1.939E-04 

4.998E-05 

-5.998E+03 

18 

4.217E-03 

-8.072E-05 

2.018E-05 

-2.422E+03 

19 

1.748E-03 

7.519E-05 

-1.880E-05 

2.256E+03 

20 

-2.664E-04 

1.998E-04 

-4.995E-05 

5.994E+03 

21 

-l.lOOE-03 

1.333E-04 

-4.982E-05 

5.978E+03 

22 

-7.572E-04 

6.927E-05 

-1.732E-05 

2.078E+03 

23 

5.487E-06 

-7.986E-05 

1.397E-05 

-2.39BE+03 

24 

3.434E-04 

-S.ESEE-OE 

2.424E-05 

-2.909E+03 

25 

1.365E-04 

-3.579E-06 

8.949E-07 

-1.074E+02 

26 

-1.137E-04 

5.467E-05 

-1.367E-05 

1.B40E+03 

27 

-7.654E-05 , 

1.603E-05 

-4.007E-06 

4.808E+02 

28 

5.425E-05 

-‘3.428E-05 

8.571E-06 

-1.029E+03 

29 

2.033E-05 

-2.156E-06 

5.389E-07 

-B.4B7E+01 

30 

-3.311E-05 

2.110E-05 

-5.275E-06 

6.330E+02 

31 

-3.5S2E-15 

-3.037E-05 

7.743E-06 

-9.291E+02 
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Fig. A-3 Displacement profiles at various times after impact of 
the steel beam. 
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Listing of Program 


PROGRAM MAIM (INPUT. 0UTPUT.PL0T»TAPE5=IMPUT»TAPEB=0UTPUTfTAPEll.TA A 
1PE8) A 

A 

CONTROL MAIN PROGRAM A 

A 

COMMON /CONTR/ TITLE(?),NP.NE,NB,NDF,NCN.NLD.NMAT,NSZF,LI,NT4,NI1IN A 
l.MATP.NPROB A 

COMMON /TIME/ T.DT.DDT.TAU.KCON.KCNT A 

COMMON /DISP/ 01, 02, 03. 010,020.030 A 

COMMON /DIME/ TB, WB, PB, NQ.Dll A 

COMMON /SPHERE/ STF,R,CABU(10),QKONST(10) A 

COMMON /PLASTIC/ DISPEM,NDISPEM,FORSPM,DISPM,OP A 

COMMON CORD(100,2),NOP(200,4),IMAT(200),ORTC25,5),NBC(25),NFIX(25) A 
1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 A 
200, 15),SK(200, 15), ISPC200,15),SMPEM(200,15),ESTIF(12,12),ENASS(12, A 
312),NFIXK(25) A 

COMMON /COMP/ 0BR(3.3,25),ABD(B,S),TH(25),ZK(25),MLAVER A 

COMMON /PLOT/ NN,TT(25),FF(25),W(25),U(25) A 

A 

INITIALIZE TAPE NO. A 

AND NUMBER OF CORNER NODE MAX. A 

A 

NT4=11 A 

NCN=2 A 

NN=1 A 

A 

PROBLEM IDENTIFICATION A 

A 

CALL PLOTS A 

101 READ (5,108) NPROB, (TITLE(I), 1=1,7) A 

IF (NPROB. EQ.O) GO TO 105 A 

DO 102 KG=1,200 A 

R10(KG)=0. A 

R20(KG)=0. A 

R30(KG)=0. A 

102 R3(KG)=0. A 

A 

READ INPUT GEOMETRY AND PROPERTIES A 

A 

CALL GDATA A 

NDISPEM=0 A 

T=0. A 

TAU=2. A 

KCON=0 A 

DDT=DT*DT A 

A 

LOOP ON NO OF PROBLEMS A 

A 

REWIND NT4 A 

NSZF=NP*NDF A 

CALL FORMK A 

CALL FORMM A 

DO 103 LI=1,NLD A 

KCNT=1 A 

A 

READ LOADS A 

A 

CALL LOAD A 

A 

FORM THEN SOLUE SIMULTANEOUS EQUATIONS A 

A 

CALL HMTO A 

CALL SOLUE A 

CALL INTEGTN A 

A 

ITERATION 2 A 

A 

KCNT=2 A 

CALL LOAD A 

CALL HMTQ A 


2 

3 

4 

5 
B 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IB 

17 

18 

19 

20 
21 
22 

23 

24 

25 
EG 

27 

28 

29 

30 

31 

32 

33 

34 

35 
3B 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 
5B 

57 

58 

59 
BO 
B1 
B2 
B3 
B4 
B5 
BB 
B7 
G8 
G9 

70 

71 
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CALL SDLUE ft 72 

CALL INTEGTN A 73 

T=T+DT ft 74 

IF (T.GT.IOO.E-G) GO TO 104 A 75 

IF (LI. EQ. 10,000) GO TO 104 ^ A 7G 

103 CONTINUE A 77 

104 WRITE (G.IOG) A 78 

WRITE (G.107) (CTT(I),FF(I),Wa)fU(I)).I=l,NN) A 79 

WRITE (8,107) ((TT(I),FF(I),W(I),U(I)),I=1,NN) A 80 

CALL FACTOR (0.8) A 81 

CALL PLOT (0.0, 2. 0,3) A 82 

CALL SCALE (TT,G.0,21, 1) A 83 

CALL SCALE (FF,9.0,21, 1) A 84 

CALL SCALES (9.0,W,21, 1,U,21, 1) A 85 

A 8B 

W(22)=U(22)=TT(22)=FF(22)=0.0 A 87 

TT(23)=20. A 88 

FF(23)=20. A 89 

W(23)=U(23)=0.001 A 90 

A 91 

CALL AXIS (0.0,0.0,10HTIME( SEC),-10,G.0,0.0,TT(22),TT(23)>0) A 92 

CALL AXIS (0.0,0.0,9HFORCE(LB),9,9.0,90.0,FF(22),FF(23),-1) A 93 

CALL AXIS (B.0,0.0,14HDISP(0.001 IN), 14,9. 0,90.0, W(22),W(23),-1) A 94 

CALL LINE (TT,FF,21, 1,1,3) A 95 

CALL DSHLINE(TT,W, 21, 0.1, 0.1,1) A 9B 

CALL DSHLINE (TT, U,21, 0. 05, 0. 05, 1) A 97 

CALL PLOT (G. 0,9. 0,3) A 98 

CALL PLOT (0.0, 9. 0,2) A 99 

CALL SYMBOL ( 1 . 0, 9.3, 0. 1, TITLE, 0. 0, 70) A 100 

CALL SYMBOL (3.5, 8.5, 0. 1, 17HBALL DIA,=l/2 IN. ,0.0, 17) A 101 

GO TO 101 A 102 

105 CALL PLOT (0,0,999) A 103 

STOP A 104 

A 105 

lOG FORMAT ( IHl , 4X, lOHTIME(MSEC) , GX, 9HF0RCE(LB) , EX, 13HBALL DISP( IN), EX A lOG 

1,13HBEAM DISP(IN)) A 107 

107 FORMAT (4E15.3) A 108 

108 FORMAT (I2,7A10) A 109 

A 110 

END A 111 

SUBROUTINE GDATA B 2 

COMMON /CONTR/ TITLE(7),NP,NE,NB,NDF,NCN,NLD,NMAT,NSZF,LI,NT4,NDIN B 3 

l,MATP,NPROB B 4 

COMMON /TIME/ T,DT, DDT, TAU,KCON,KCNT B 5 

COMMON /DISP/ 01,02,03,010,020,030 B G 

COMMON /DIMB/ TB, WB, PB, NO, D1 1 B 7 

COMMON /SPHERE/ STF,R,CABU(10),OKONST(10) B 8 

COMMON /PLASTIC/ DISPEM, NDISPEM, FORSPM, DISPM, OR B 9 

COMMON CDRD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25) B 10 

1,R1(200),R2(200),R3(EOO),R10(200),R20(200),R30(200),FORS(EOO),SM(2 B 11 

200,15),SK(200,15),ISP(200,15),SMPEM(E00,15),ESTIF(12,12),EMASS(12, B 12 

312),NFIXK(25) B 13 

COMMON /COMP/ QBR(3,3,25),ABD(G,G),TH(E5),ZK(25),MLAYER B 14 

COMMON /PLOT/ NN, TT(25) , FF(25) , W(E5) , U(25) B 15 

B IB 

READ AND PRINT TITLE AND CONTROL B 17 

B 18 

WRITE (G,11G) NPROB, (TITLE(I), 1=1,7) B 19 

WRITE (8,11B) NPROB, (TITLE(I), 1=1,7) B 20 

READ (5,106) NP,NE,NB,NTM,NMAT,NDIN,MATP,NDC,I1 B 21 

NDF=3 B EE 

NLD=NDIN»NTM B 23 

FLD=FLDAT(NLD)/10. B 24 

FDIN=FLDAT(NDIN)/10. B 25 

READ (5,113) TB,WB,R,N0,02,DT B EG 

WRITE (G,107) NP,NE,NB,NLD,NDF,NMAT B 27 

NLD=NLD+1 B 28 

B 29 

READ AND PRINT MATERIAL DATA B 30 

SPHERE DATA! L=NMAT (LAST MAT. CARD) B 31 
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C B 32 

READ (5.112) (L,C0RT(L,I).;I=1.5).M=l.NmT) B 33 

PB=0RT(NMAT,5) B 34 

WRITE (B.122) TB.WB.PB.R.NQ.QE.DT B 35 

NQ=(NQ-1)«3+1 B 3B 

WRITE (B. 121) B 37 

WRITE (B.115) (N. (DRTCH. I).I=1.5).M=l,mftT) B 38 

B 39 

READ INDENTATION DATA B 40 

B 41 

READ (5.111) STF.DISPEM.QP B 42 

WRITE (B.123) DISPEM B 43 

IF (DISPEN.NE.0.0) WRITE (B.124) QP B 44 

B 45 

READ NODAL POINT DATA B 4B 

AND B 47 

READ ELEMENT DATA B 48 

B 49 

DO 102 1=1. NDC B 50 

READ (5.114) ND1.ND2.X1.X2.IMT B 51 

EL=(X2-X1)/FL0AT(ND2-ND1) B 52 

CORD(NDl. 1)=X1 B 53 

C0RD(ND2. 1)=X2 B 54 

CORD(ND2.2)=0.() B 55 

C0RD(ND1.2)=C0RD(ND2.2) B 5B 

NDD=ND2-1 B 57 

DO 101 J=ND1.NDD B 58 

CORD(J+l.l)=CORD(J.l)+EL B 59 

C0RD(J+1.2)=0.0 B BO 

N0P(J,1)=J B B1 

NDP(J.2)=J+1 B B2 

NDP(J.4)=0 B B3 

N0P(J.3)=N0P(J.4) B 64 

IMAT(J)=IMT B B5 

101 CONTINUE B 6B 

102 CONTINUE B 67 

B B8 

READ BOUNDARY DATA B 69 

B 70 

READ (5.110) (NBC(I).NFIXCI).I=1.NB) B 71 

IF (MATP.EQ.l) CALL CliPD B 72 

B 73 

ISOTROPIC MATP=0.0 B 74 

COMPOSITE MATP=1.0 B 75 

B 7S 

IF (Il.NE.O) GO TO 103 B 77 

B 78 

PRINT INPUT DATA B 79 

B 80 

WRITE (6.117) B 81 

WRITE (B.108) (N.(C0RD(N.M).M=1.2).N=1.NP) B 82 

WRITE (B.118) B 83 

WRITE (B.109) (N. (N0P(N.M).M=1.4).IMAT(N).N=1.NE) B 84 

WRITE (6.119) B 85 

WRITE (6.110) (NBC(I).NFIX(I). I=1.NB) B 86 

WRITE (6.120) FDIN.FLD B 87 

103 CONTINUE B 88 

DO 104 IJ=1.200 B 89 

R10(IJ)=0. B 90 

R20(IJ)=0. B 91 

R3O(IJ)=0. B 92 

104 FORS(IJ)=0. B 93 

DO 105 IJ=1.25 B 94 

105 NFIXK(IJ)=NFIX(IJ) B 95 

RETURN B 96 

C B 97 

106 FORMAT (915) B 98 

107 FORMAT (13HONODAL POINTS. 9X. I5/1X.8HELEMENTS, 13X. I5/1X. 19HB0UNDARY B 99 

1 CONDITIONS. EX. I5/1X.12H0UTPUT LIMIT. lOX. I5/1X. 18HDEGREES OF FREED B 100 

20M.3X.I5/1X.9HMATERIALS.12X.I5) B 101 



108 FORMAT (110,2F10.3) 

109 FORMAT (BIS) 

110 FORMAT (215) 

111 FORMAT (E10.3.2F10.0) 

112 FORMAT (I5»5F10.0) 

113 FORMAT OFIO.O/IS.RFIO.O) 

114 FORMAT (215. EFIO. 0.15) 

115 FORMAT (I5.7X.3(F10.1.4X)1F5.3,7X.F8.B//) 

IIG FORMAT (1H1.I2.7A10) 

117 FORMAT (14H0 MODAL POINTS/'ITX. IHX. lOX. IHY) 

118 FORMAT (lOHO ELEMENTS/9X.1HI.4X.1HJ.4X.1HK.8X.3HMAT) 

119 FORMAT (21H0 BDUMDARY CONDITIONS) 

120 FORMAT (IBHOPRINTING 5CHEME/5X.22H1. REPORT OUTPUT EUERY.FB.2.2X.4 
1HMSEC/5X.22H2. TERMINATE OUTPUT AT.FB.2.2X.4HMSEC) 

121 FORMAT (1H0.20H MATERIAL PR0PERTIES/1X.8HMAT. N0..7X.2HE1.12X.2HE2 
1 . 1 IX. 3HG12. 1 OX, 3HU12, 1 OX, 3HRH0/) 

122 FORMAT (15H0BEAM THICKNESS, 1 IX. FB.3/1X, lOHBEAM WIDTH, 15X.FB.3/1X, 1 
14HSPHERE DENSITY, 12X.F8.G/1X, 13HSPHERE RADIUS, 12X,FG.3//1X, IIHIMPA 
2CT NODE, 14X, I2/-1X, 15HIMPACT UELOCITY, 10X.F6.1/1X, 39HINTEGRATI0N T 
3IME INCREMENT ( X E-OB SEC), El 0.3) 

123 FORMAT (//,1X, 21HPERMANENT DEFORMATION, 9X.F8.B) 

124 FORMAT (/, IX, 15HUNL0ADING POWER, 15X,FB.3) 


END 

SUBROUTINE ESTIFM (N) 

REAL IB, LB 

COMMON /CONTR/ TITLE (7) , NPi NE, NB, NDF', NCN, NUD, NMAT, NS2F, LI ii.NT4, NDIN 
l.MATP.NPROB 

COMMON /'TIME/' T.DT.DDT.TAU.KCON.KCNT 
COMMON /DISP/' 01,02,03,010,020,030 
COMMON /'DIME/' TB,WB,PB,NQ,:D11 

COMMON CORD(100,2),NOP(200,4),IMAT(2D0),ORT(25,6),NBC(25),NFIX(25) 
1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(I200),FORS(200),SM(2 
200,15), SK(200, 15) , ISP(200, 15) , SMPEM(200, 15) , ESTIF( 12, 12) , EMASS( 12, 
312),NFIXK(25) 

COMMON /'COMP/' 0BR(3,3,25)iABD(S,G),TH(25),ZK(25),MLAYER 
IB=WB*TB*»3/12. 

LB=C0RD(N+1, 1 )-CORD(N, 1) 

SOLB=LB*LB 

TPLB=LB*LB*LB 

1MN=IMAT(N) 

PARA1=0RT( IMN, 1 )*IB/70 . 

IF (MATP.EO.l) PARAl*ABD(4,4)/70. 

ESTIF ( 1 , 1 )=1200 . /TPLB^tPARAl 
ESTIF ( 1 , 2) =B00 .i/SOLB*PARAl' 

ESTI F ( 1 , 3)=30 . /LB*PARA1 
ESTIF (1,4) =-1200. /TPLB»PARA1 
ESTIF( 1 , 5)=B00 .i/'SQLB*PARAl 
ESTIF ( 1 . B ) =-30 ../LB*PARA1 
ESTIF(2,1)=ESTIF(1,2) 

EST I F ( 2 , 2 ) =384 .1/'LB*PARA1 
ESTIF(2,3)=22.«*PARA1 
ESTIF( 2, 4 )=-BOQ . /SQLB*PARA1 
ESTIF(2,5)=21B.i/LB*PARAl - 
ESTIF(2,B)=-8.*PARA1 
ESTIF(3,1)=ESTIF(1,3) 

ESTIF(3,2)=ESTIF(2,3) 

ESTIF(3,3)=B.»LB»PARA1 
ESTIF(3,4)=-30.;/LB*PARAl 
ESTIF ( 3 , 5 ) =8 . *PARA1 
ESTIF(3,B)=LB*PARA1 
ESTIF(4,1)=ESTIF(1,4) 

ESTIF(4,2)=ESTIF(2,4) 

ESTIF(4,3)=ESTIF(3,4) 

ESTIF(4,4)=1200./TPLB*PARA1 
ESTIF( 4, 5 )=-B00 . /'S0LB*PARA1 
ESTIF(4, B )=30 . /LB*PARA1 
ESTIF(5,l)=ESTrF(1.5) 

ESTIF(5,2)=ESTIF(2,5) 

ESTIF(5,3)=ESTIF(3,5) 


B 102 
B 103 
B 104 
B 105 
B 106 
B 107 
B 108 
B 109 
B 110 
B 111 
B 112 
B 113 
B 114 
B 115 
B 116 
B 117 
B 118 
B 119 
B 120 
B 121 
B 122 
B 123 
B 124 
B 125 
C 2 
C 3 
C 4 
C 5 
C *6 
C 7 
C 8 
C 9 
C 10 
C 11 
C 12 
C 13 
C 14 
C 15 
C 16 
C 17 
C 18 
C 19 
C 20 
C 21 
C 22 
C 23 
C 24 
C 25 
C 26 
C 27 
C 28 
C 29 
C 30 
C 31 
C 32 
C 33 
C 34 
C 35 
C 36 
C 37 
C 38 
C 39 
C 40 
C 41 
C 42 
C 43 
C 44 
C 45 
C 4B 
C 47 
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ESTIF(5.4)=ESTIF(4»5) C 

ESTIF(5,5)=384./LB«PftRftl C 

ESTIF(5,G)=-22.«PARA1 C 

ESTIF(B, l)=ESTI.Fa,G) C 

ESTIF(B,2)=ESTIF(2»B) C 

ESTIF(B,3)=ESTIF(3.6) C 

ESTIF(6»4)=ESTIF(4»6) C 

ESTIF(B,5)=ESTIF(5!.B) C 

ESTIF(B,B)=B.*LB«PARA1 C 

IF (N.NE.l) GO TO 101 C 

WRITE (B.103) C 

WRITE (Ga02) ((ESTIFa,J)»J=l,B).I=l.B) C 

101 CONTINUE C 

RETURN C 

c c 

102 FORMAT (1X,BE11.3) C 

103 FORMAT (1H0>38H TYPICAL STIFNESS MATRIX OF AN ELEMENT) C 

C C 

END C 

SUBROUTINE EMASSM (N) D 

REAL LB D 

COMMON /CONTR/ TITLE(?),NP,NE,NB»NDF.NCN,NLD,NMAT.NSZF,LI.NT4»NDIN D 
l.MATP.NPRDB D 

COMMON /TIME/ TpDTsDDTaAU.KCON.KCNT D 

COMMON /DISP/ Q1.Q2, 03, 010.020,030 D 

COMMON /DIME/ TB,WB,PB,NQ,D11 D 

COMMON CORD(100,2),NOP(200,4),IMATC200),ORT(25,5),NBC(25),NFIX(25) D 
1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 D 
200,15),SK(200,15),ISP(200,15),SMPEM(200.15),ESTIF(12,12),EMASS(12, D 
312),NFIXKC25) D 

COMMON /COMP/ QBR(3,3,25),ABD(G,B),TH(25),ZK(25),MLAYER D 

LB=C0RD(N+1, 1)“C0RD(N, 1) D 

AB=TB»WB D 

SQLB=LB*fLB D 

TPLB=LB-»LBi^LB D 

QDLB=LB-»LB^^LB«LB D 

IMN=IMAT(N) D 

PARA2=ORTaMN,5)*AB»LB/55440. D 

EMASSCl, 1)=21720.«PARA2 D 

EMASSC 1,2) =3732 aLB*PARA2 D 

EMASS(1,3)=281.«SQLB«PARA2 D 

EMASS(1,4)=B000.»PARA2 D 

EMASS(l,5)=-1812.»LB»fPARA2 D 

EMASS(1,B)=181.«SQLB»PARA2 D 

EMASS(2, 1)=EMASS(1,2) D 

EMASS(2,2)=832.*SQLB«PARA2 D 

EMASS(2,3)=B9aTPLB*PARA2 D 

EMASS(2,4) = 1812.»tLB^fPARA2 D 

EMASS(2,5)=-532.«SQLB«PARA2 D 

EMASS(2,B)=52.«-TPLB*PARA2 D 

EMASS(3,1)=EMASS(1,3) D 

EMASS(3,2)=EMASS(2,3) D 

EMASS(3,3)=G.«QDLB»PARA2 D 

EMASS(3,4)=181.»SQLB*PARA2 D 

EMASS(3,5)=-52.*TPLB»PARA2 D 

EMASS(3,B)=5.*QDLB»PARA2 D 

EMASS(4,1)=EMASS(1,4) D 

EMASS(4,2)=EMASS(2,4) D 

EMASS(4,3)=EI1ASS(3,4) D 

EMASS(4,4)=21720.»PARA2 D 

EMASS(4,5)=-3732.»LB^fPARA2 D 

EMASS(4,G)=281.*SQLB^tPARA2 D 

EMASS(5, 1)-EMASS(1,5) D 

EMASS(5,2)=EMASS(2,5) D 

EMASS(5,3)=EMASS(3,5) D 

EMASS(5,4)=EMASSC4,5) D 

EMASS(5,5)=832.*S0LB»PARAE D 

EMASS(5,B)=:-G9.»TPLB*PARA2 D 

EMASS(B,1)=EMASS(1,B) D 

EMASS(B,2)-=EMA5S(2,B) D 
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EMASS(G,3)=EMflSS(3,G) B 53 

EMASS(G,4)=EM(=lSS(4,G) D 54 

EMASS(G,5)=EMASS(5,G) D 55 

EmSS(G,G)=6.*QDLB*PARA2 D 5G 

IF (N.NE.l) GO TO 101 D 57 

WRITE (G,103) D 58 

WRITE (Gp102) ((EMASS(I,J),J=1,S),I=1,B) D 59 

101 CONTINUE D BO 

RETURN B B1 

C B G2 

102 FORMAT (lXiGEll.3) B G3 

103 FORMAT (1H0,34H TYPICAL MASS MATRIX OF AN ELEMENT) D 64 

C D G5 

ENB B BB 

SUBROUTINE FORMM E 2 

E 3 

FORMS MASS MATRIX E 4 

IN UPPER TRIANGULAR FORM E 5 

E B 

COMMON /CONTR/ TITLE(7),NP,NE»NB,NBFpNCN,MLBpNMAT,NSZF.LI,.NT4»NBIN E 7 

IpMATP.NPROB E 8 

COMMON /TIME/ T, BT. BBT, TAU» KCONp KCNT E 9 

COMMON /BISP/ 01, 02, 03, 010,020,030 E 10 

COMMON /BIMB/ TB, WB, PB,N0,B11 E 11 

COMMON CORB(100,2),NOP(200,4),IMAT(200),DRTC25,5),NBC(25),NFIX(25) E 12 

1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 E 13 

200,15),SK(200,15),ISP(200,15),SMPEM(200,15),ESTIF(12,12),EMASSa2, E 14 

312),NFIXK(25) E 15 

COMMON /COMP/ QBRi:3,3,25),ABB(B,6),TH(25),ZK(25),MLAYER E 16 

E 17 

SET BANBMAX ANB NO. OF EOUATIONS E 18 

E 19 

NBANB=9 E 20 

E 21 

ZERO MASS MATRIX E 22 

E 23 

BO 101 N=1,NSZF E 24 

BO 101 M=1,NBANB E 25 

101 SM(N,M)=0. E 2B 

E 27 

SCAN ELEMENTS E 28 

E 29 

BO lOB N=1,NE E 30 

CALL EMASSM (N) E 31 

E 32 

RETURNS EMASS AS MASS MATRIX E 33 

E 34 

STORE EMASS IN SM E 35 

E 3G 

FIRST ROWS E 37 

E 38 

BO 105 JJ=1,NCN E 39 

NROWB=(NOP(N, JJ)-1)*NBF E 40 

BO 105 J=1,NBF E 41 

NR0WB=NR0WB+1 E 42 

I=(JJ-1)«NBF+J E 43 

E 44 

THEN COLUMNS E 45 

E 4G 

DO 104 KK=1,NCN E 47 

NC0LB=(N0P(N,KK)-1)»NBF E 48 

BO 103 K=1,NBF E 49 

L=(KK-1)*NBF+K E 50 

NC0L=NC0LB+K+1-NR0WB E 51 

E 52 

SKIP STORING IF BELOW BANB E 53 

E 54 

IF (NCOL) 103,103,102 E 55 

102 SM(NROWB,NCOL)=SM(NROHB,NCOL)+EMASS(I,L) E 5G 

103 CONTINUE E 57 
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104 CONTIMUE 

105 CONTINUE 
lOG CONTINUE 

C 

C INSERT BOUNDARY CONDITIONS 

C 

DO 112 N=1.NB 
NX=10*-»(NDF-1) 

I=NBC(N) 

NROUB=CI-l)*NDF 

EXAMINE EACH DEGREE OF FREEDOM 

DO 111 M=1,NDF 
NROWB=NROWB+l 
ICON=NFIX(N)/NX 
IF (ICON) 110i.ll0»107 
SMCNROUB, 1)=1. 

DO 109 J=2,NBAND 
SMCNROWB, J)=0. 

NR=NROWB+l-J 
IF (NR) 109,109*108 
SM(NR, J)=0. 

CONTINUE 

NFIX(N)=NFIX(N)-NX*ICON 
NX=NX/10 
CONTINUE 

112 CONTINUE 

DO 115 N=1,NSZF 
K=0 

DO 114 M=1,NBAND 
MP=M-K 

IF (ISP(N,M).LT.ISP(N,D) GO TO 113 
SM(N,MP)=SM(N,MP)+(DDT/G.)*SK(N,M) 

GO TO 114 

113 K=K+1 

114 CONTINUE 

115 CONTINUE 

DO lie I=1,NSZF 
DO IIB J-1,NBAND 
IIG SMPEM(I,J)=SM(I,J) 

C 

C WRITE(G,1) ((SM(I, J), J=1,NBAND),I=1,NSZF) 

C 1 FQRMAT(2X,9E10.3) 

C 

RETURN 

C 

END 

SUBROUTINE FDRMK 

COMMON /CDNTR/ T I TLE C 7 ) , NP , NE, NB , NDF , NCN, NLD, NMAT , NSZF , LI ,.NT4, NDIN 
l.MATP.NPROB 

COMMON /TIME/ T, DT, DDT, TAU, KCON, KCNT 
COMMON /DISP/ 01,02, 03, 010, 020,030 
COMMON /DIME/ TB, WB, PB, NO, Dll 

COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(E5) 
1,R1(200),R2(EOO),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 
S00,15),SK(800,15),ISP(200,15),SMPEM(200,15),ESTIF(1E,1S),EMASS(1E, 
31E),NFIXK(E5) 

COMMON /COMP/ QBR(3,3,S5),ABD(G,G),TH(E5),ZK(E5),MLAYER 
C 

C SET MAX. NO. OF TERMS 

C 

NMAX=9 

N0FF=9 

C 

C ZERO ARRAYS 

C 

DO 103 N=1,NSZF 
DO 101 M=1,NMAX 
101 SK(N,M)=0. 
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111 


E 58 

E 59 

E GO 

E 61 

E BE 

E G3 

E B4 

E G5 

E GG 

E B7 

E B8 

E G9 

E 70 

E 71 

E 7E 

E 73 

E 74 

E 75 

E 7G 

E 77 

E 78 

E 79 

E 80 

E 81 

E 8S 

E 83 

E 84 

E 85 

E 8B 

E 87 

E 88 

E 89 

E 90 

E 91 

E 92 

E 93 

E 94 

E 95 

E 96 

E 97 

E 98 

E 99 

E 100 

E 101 

E 102 

E 103 

E 104 

E 105 

F 2 

F 3 

F 4 

F 5 

F B 

F' 7 

F 8 

F 9 

F 10 

F 11 

F 12 

F 13 

F 14 

F 15 

F IB 

F 17 

F 18 

F 19 

F 20 

F 21 

F 22 

F 23 



91 


DO 102 M=2,N0FF F 24 

102 ISP(N.M)=0 F 25 

103 ISP(N,l)=rS F 2G 

C F 27 

C SCAN ELEMENTS F 28 

C F 29 

DO 110 N=1.NE F 30 

CALL ESTIFM (N) F 31 

C F 32 

C RETURNS ESTIF AS STIFFNESS MATRIX F 33 

C F 34 

C STORE ESTIF IN SK NITH A TERM IN ISP AS A POINTER F 35 

C F 3B 

C F 37 

C FIRST THE ROWS F 38 

C F 39 

1=0 F 40 

DO 109 JJ=1.NCN F 41 

NROUB=(NOP(N, JJ)-1)«NDF F 42 

DO 109 J=1,NDF F 43 

NRQWB=NROWB+l F 44 

1=1+1 F 45 

C F 48 

C THEN COLUMNS OF ESTIF F 47 

C F 48 

11=0 F 49 

DO 108 KK=1,NCN F 50 

NCOLB=(NOP(N.KK)~l)*NDF F 51 

DO 108 K=1.NDF F 52 

NCOLB=NCOLB+l F 53 

11=11+1 F 54 

C F 55 

C SEARCH ISP FOR COLUMN NO. F 58 

C F 57 

DO 105- M=l.NOFF F 58 

IF (ISP(NROWB,M)-NCOLB) 104,107»104 F 59 

104 IF (ISPCNROUB.M)) 10G»10B»105 F 80 

105 CONTINUE F 81 

C F 82 

C FOUND A BLANK NOW STORE NCOLB F 83 

C F 84 

108 ISP(NROWB.M)=NCOLB F 85 

C F 88 

C NOW STORE ESTIF F 67 

C F 88 

107 SK(NROWB,M)=ESTIFa,II)+SK(NROWB,M) F 89 

C F 70 

C END LOOP ON COLUMNS F 71 

C F 72 

108 CONTINUE F 73 

C F 74 

C END LOOP ON ROWS F 75 

C F 78 

109 CONTINUE F 77 

C F 78 

C END LOOP ON ELEMENTS F 79 

C F 80 

110 CONTINUE F 81 

C F 82 

C INSERT BOUNDARY CONDITIONS F 83 

C F 84 

)DO 114 N=1,NB F 85 

NX=10»»(NDF--1) F 88 

I=NBC(N) F 87 

NROWB=(I-l)^fNDF F 88 

C F 89 

C EXAMINE EACH DEGREE OF FREEDOM F 90 

C F 91 

DO 113 M=1,NDF F 92 

NROWB=NROWB+l F 93 
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ICOM=IHFI><K(M)/M><! F 

IF (ICON) 112,112.111 F 

F 

STORE ZERO ON DIAGONAL F 

F 

111 SK(NR0WB, 1)=0.0 F 

NFIXK(N)=NFIXK(N)-NX*ICON F 

112 NX=NX/10 F 

113 CONTINUE F 

114 CONTINUE F 

RETURN F 

F 

END F 

SUBROUTINE LOAD G 

COMMON /CONTR/ TITLE(7),NP,NE,NB»NDF.NCN.NLD.NNAT.NS2F.LI,NT4.NDIN G 
l.MATP.NPROB G 

COMMON /TIME/ T, DT, DDT, TAU, KCON. KCNT G 

COMMON /DISP/ Q1.Q2, 03, QIO, 020,030 G 

COMMON /DIME/ TB, WB, PB, NQ,D11 G 

COMMON /SPHERE/ STF,R,CABU(10),OKONST(10) G 

COMMON /PLASTIC/ DISPEM,NDISPEM,FORSPM,DISPM,QP G 

COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25) G 
1,R1(200),R2(200),R3(200),R10C200),R20(200),R30(200),FDRS(200),SMC2 G 

200, 15),SK(200, 15), ISPC200, 15),SMPEM(200, 15), ESTIFC 12, 12), EMASSC 12, G 

312),NFIXK(2S) G 

COMMON /COMP/ QBR(3,3,25),ABD(B,B),TH(25),ZK(E5),MLAYER G 

IF (STF.NE.0.0) GO TO 101 G 

STFI=(4./3.)*SQRT(R)/((1.-0RT(NMAT,4)»*2)/DRT(NMAT, l) + a.-0RTa,4) G 
I«<f2)/0RT(1, D) G 

STF A= ( 4 . /3 . ) «SQRT ( R ) / ( ( 1 . -ORT ( NMAT, 4 ) ^f*2 ) /ORT ( NMAT , 1 ) + 1 . /ORT (1,2)) G 
STF=STFI G 

IF (MATP.EQ.l) STF=STFA G 

101 PAI=4,«ATAN(1. ) G 

BALLM=(4./3. )*PAI»(Rit^«3)<«PB G 

G 

SIMPLY SUPPORTED SYMMETRY CST1=0.5 G 

CLAMPED CANTILEUER CST1=1. G 

G 

G 

CST1=1.0 G 

G 

IF(NBCd) .EQ. 1) CST1=0.5 G 

G 

Q1=ACCEL. OF^ THE MASS G 

Q2=UEL0. OF THE MASS G 

Q3=DISP. OF THE MASS G 

G 

IF (LI.GT.l.AND.KCNT.EO.l) GO TO 102 G 

IF (LI.GT,1.AND.KCNT.EQ.2) GO TO 103' G 

0i=0. . G 

03=0, G 

GO TO 112 G 

102 010=01 G 

020=02 G 

030=03 G 

Q3=030+DT*tiEO+0,5*DDT«010 G 

R3(NQ)=R30(NQ)+DT*P20(HQ)+0.5*DDT*R10(NQ) G 

DIFD0=Q30-R30(NQ) G 

DIFDISP=Q3-R3(NQ) G 

G 

WRITE(B,400) DIFDO.DIFDISP G 

G 

IF (DIFDISP) 110,104,104 G 

103 a3=Q30+DT»Q20+DDT*Q10/3.+DDT*Ql/G, G 

DIFDO=Q30-R3O(NO) G 

D1FD1SP=Q3-R3(NQ) G 

G 

WR1TE(G,400) DIFDO.DIFDISP G 

400 F0RMAT(/,5X,?iDlFD0=?^,E15.3,5X,?«DIFDISP=?!,E15.3) G 

G 
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17 
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19 
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21 
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IF (DIFDISP.lt. 0) GO TO 110 G 59 

104 IF (DISPEM.EQ.0.0) GO TO 105 G BO 

IF ((DIFDISP.LT.DIFDD).AMD.(NDISPEM.EQ.O)) GO TO 107 G B1 

IF ((DIFDISP.LT.DIFDO).AND.(NDISPEM.GT.O)) GO TO 108 G G2 

105 DO lOB J=1,NSZF , G B3 

lOG FORS(J)=0. G 64 

F0RS(NQ)=STF»(DIFDISP)»»1.5*CST1 G 65 

Q1=-F0RS(NQ)/BALLM/CST1 G GG 

IF (KCNT.EQ.l) GO TO 113 G G7 

IF (KCNT.EQ.2) GO TO 109 G 68 

107 NDISPEt1=l G G9 

FORSPM=FDRS(ISQ) G 70 

DISPM=DIFDO G 71 

WRITE (S.114) DISPEM.DISPM.DIFDISP.FORSPM G 72 

IF ((DIFDISP.LT.DISPEM).OR.(DISPM.LE.DISPEM)) GO TO 111 G 73 

108 FORS(MQ)=FORSPM*((DIFDISP-DISPEM)/(DISPM-DISPEM))**QP*CSTl G 74 

Ql=-FORS(NQ)/BALLM/CSTl G 75 

IF (KCMT.EQ.l) GO TO 113 G 7G 

109 Q2=Q20+0.5*DT»(D10+0.5*DT*Q1 G 77 

Q3=Q30+DT*Q20+DDT*Q10/3.+DDT*Ql/G. G 78 

GO TO 113 G 79 

110 FORS(MQ)=0. G 80 

Q1=0. G 81 

GO TO 109 G 82 

111 LI=10000 G 83 

GO TO 113 G 84 

112 FORS(MQ)=0. G 85 

113 RETURM G 8B 

C G 87 

114 FORMAT C///,5Xf 7HDISPEM=,E10.3.5X. BHDISPM=. E10.3. 5X. 8HDIFDIS G 88 

1P=.E10.3.5X. 7HF0RSPM=,E10.3) G 89 

C G 90 

END G 91 

SUBROUTIME HMT(3 H 2 

H 3 

SUBROUTIME F'OR FIMDIMG (F)-(K)(U) H 4 

H 5 

COMMON /COMTR/ TITLE(7).MP.NE,NB,NDr,NCN,NLD»NMAT,NSZF,LI*;NT4,NDIN H G 

l,MATP.NPROB H 7 

COMMON /TIME/ T. DT, DDT, TAU, KCON. KCNT H 8 

COMMON /DISP/ CJ1,Q2,Q3,Q10,Q20,Q30 H 9 

COMMON /DIME/ TB,WB,PB,NQ,‘D11 H 10 

COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(E5) H 11 

1,R1(200),R2(200),R3(200),R10(200),R20(200),R30C200),FORS(200),SM(2 H 12 

200,15),SK(200,15),ISP(200,15),SMPEM(200,15),ESTIF(12,12),EMASS(12, H 13 

312),NFIXK(25) H 14 

COMMON /COMP/ QBR(3,3,25),ABD(G,G),TH(25),2K(25),MLAYER H 15 

NT=9 H IB 

DO 101 IJ=1,NSZF H 17 

R1(IJ)=0. H 18 

R2(IJ)=0. H 19 

101 R3(IJ)=0. H 20 

H 21 

H 22 

DO 105 N=1,NSZF H 23 

FX=FORS(N) H 24 

DO 102 M=1,NT H 25 

L=ISP(N,M) H EG 

102 FX=FX-SK(N,M)«(R30(L)+DT*R20(L)+(DDT/3.0*R10(U) H 27 

IF (SKCN.D) 104,103,104 H 28 

103 FX=0. H 29 

104 R1(N)=FX H 30 

105 CONTINUE H 31 

RETURN H 32 

H 33 

END H 34 

SUBROUTINE SOLUE I 2 

I 3 

SPECIFICATION STATEMENTS I 4 

I 5 
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COtIMOM /CONTR/ TITLE(7)!>NP.NE,tiB.raF»NCN»MLD»MMAT.NSZF»LI»NT4,NDIM I 
l.MATPfNPRDB I 

COMMOM /TIME/ T, DT, DDT. TAU, KCOM, KCNT I 

COMMON /DISP/ 01,02.03.010.020,030 I 

COMMON /DIME/ TB.UB.PB.NQ.Dll I 

COMMON CORDa00.2),NOP(200.4).IMAT(200),ORT(25,5),NBCC25),NFIX(25) I 
l.Rl(200),R2(200).R3(200).R10(200),R20(200),R30(200).FORS(200),SM(2 I 

200.15) ,SK(200,15).ISP(200,15).SMPEM(200,15),ESTIF(12,12),EMASS(12, I 

312),NFIXK(25) I 

COMMON /COMP/ QBR(3.3,25).ABD(G.B),TH(25).2K(25I.MLAYER I 

NBAND=9 I 

DO 101 I=1.NS2F I 

DO 101 J=1.NBAND I 

101 SM(1.J)=SMPEM(I. J) I 

I 

REDUCE MATRIX I 

I 

DO 106 N=1.NS2F I 

I=N I 

DO 105 L=2.NBAND I 

1=1+1 I 

IF (SM(N.L)) 102,105.102 I 

102 C=SM(N.L)/SM(N. 1) I 

J=0 I 

DO 104 K=L,NBAND I 

J=J+1 I 

IF (SM(N.K)) 103/104,103 I 

103 SM(I. J)=SM(I, J)-C^fSM(N.K) I 

104 CONTINUE , I 

SM(N.L)=C I 

I 

AND LOAD UECTOR I 

FOR EACH EQUATION I 

I 

R1(I)=R1(I)-C«R1(N) I 

105 CONTINUE I 

lOG R1CN)=R1(N)/SM(N,1) I 

I 

BACK-SUBSTITUTION I 

I 

N=NSZF I 

107 N=N-1 I 

IF (N) 111.111,108 I 

108 L=N I 

DO 110 K=2.NBAND I 

L=L+1 I 

IF (SM(N.K)) 109,110,109 I 

109 R1(N)=R1(N)-SM(N,K)«R1(L) I 

110 CONTINUE I 

GO TO 107 I 

111 RETURN I 

I 

END I 

SUBROUTINE INTEGTN J 

COMMON /CONTR/ TITLE(7),NP,NE,NB,NDF,NCN,NLD,NMAT,NSZF,LI,NT4,NDIN J 
l.MATP.NPROB J 

COMMON /TIME/ T. DT, DDT, TAU, KCON, KCNT J 

COMMON /DISP/ 01.02.03,010,020,030 J 

COMMON /DIME/ TB,WB,PB,NQ,D11 J 

COMMON CORDa00,8),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25) J 
l.Rl(200),R2(200).R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 J 

200.15) .SK(200.15).ISP(200, 15),SMPEM(200.15),ESTIF(12, 12).EMASS(12, J 

312),NFIXK(25) J 

COMMON /COMP/ QBR(3,3,25),ABD(G,6).TH(25),ZK(25),MLAYER J 

COMMON /PLOT/ NN,TT(25),FF(25),W(25).U(25) J 

J 

R1=ACCEL. OF BEAM J 

R2=UELD. OF BEAM J 

R3=DISPL. OF BEAM J 

J 


G 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IG 

17 

18 

19 

20 
21 
22 

23 

24 

25 
2B 

27 

28 

29 

30 

31 

32 

33 

34 

35 
3B 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 
5B 

57 

58 
2 

3 

4 

5 
B 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IG 

17 

18 
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:00 101 IJ=1»ISS2F J 19 

R2(IJ)=R20aJl+0.5*DT»R10aj)+0.5»DT«Rl(IJ) J 20 

R3aJ)=R30(;CJ)+DT»R20aJ)+(DDT/3.)*R10aJ)+CDDT/B.)*Rl(IJ) J 21 

101 COMTINUE J 29 

IF (KCHT.EQ.l) GO TO 107 J 23 

DO 102 1K=1.NP J 24 

IK4=(IK-1)»3+1 J 25 

R30aK)=R3(IK4) J 2G 

102 CONTINUE J 27 

IF ((LI/10000). EQ.l) GO TO 103 J 28 

J 29 

PRINT CONTROL J 30 

J 31 

NT0N=(LI-1)/NDIN J 32 

IF (NTON.NE.KCON) GO TO 105 J 33 

103 CONTINUE J 34 

J 35 

SINPLY SUPPORTED BEAM CST2=2. J 3G 

CANTILEUER CST2=1. J 37 

J 38 

CST2=1. J 39 

J 40 

IF(NBC(1).EQ.1)CST2=2. J 41 

J 42 

F=CST2*F0RS(NQ) J 43 

APHA=Q3-R30(NQ) J 44 

FF(NN)=F J 45 

N(NN)=Q3»1000. J 4G 

U(NN)=R30(NQ)*1000. J 47 

URITE (G,108) (TITLE(I)*I=1,7) J 48 

T1=T»1.EG J 49 

TT(NN)=T1 J 50 

NN=NN+1 J 51 

NRITE (G.109) T1.F.Q3.Q2.01,APHA J 52 

DO 104 IK=1.NP J 53 

IK3=IK*3 J 54 

ST><;X=R3(IK3)»TB/2. J 55 

SIGX=0RT(1.1)*STXX J 5G 

STYY=-0RT(1,4)»STXX J 57 

WRITE (G.llO) IK»R30(IK),STXX.STYY»SIGX J 58 

104 CONTINUE J 59 

KCON=KCON+l J GO 

105 CONTINUE J 61 

DO lOG IJ=1.NSZF J G2 

R10(IJ)=R1(IJ) J B3 

R20(IJ)=R2(IJ) J G4 

lOB R30(IJ)=R3(IJ) J 65 

107 RETURN J 6G 

C J G7 

108 F-ORMAT (1H1.7A10///) J 68 

109 F'ORMAT (10X,18HTIME ELAPSED(MSEC). 13X»F7.3/10X,9HFORCE(LB),21X,E11 J 69 

1,.3/10X»21HMASS DISPLACEMENT(IN)»9X,E11.3/10X,21HMASS UELOCITY( IN/S J 70 

2E:C),9X,E11.3/10X,20HMASS ACCEL. (IN/SEC2), lOX, Ell. 3/lOX, 15HINDENTAT J 71 

3:tON(IN), 15X,E11.3///10X,4HNODE,9X,4HDISP. 13X. 9HSTRAIN-XX. 9X, 9HSTRA J 72 

4:tN-YY,9X,9HSTRESS-XX/) J 73 

110 f-ORMAT (9X, I3,7X,E12.3f7X,E12.3»7X»E12.3,7X»E12.3) J 74 

C J 75 

END J 76 

SUBROUTINE CMPD K 2 

COMMON /CONTR/ TITLEC7) , NP, NE, MB, NDF. NON, NLD, NMAT, NS2F, LI , NT4» NDIN K 3 

l.MATP.NPROB K 4 

COMMON /TIME/ T.DT.DDT.TAU.KCON.KCNT K 5 

COMMON /DISP/ Q1,Q2.Q3,Q10.Q20»Q30 K 6 

COMMON /DIMB/ TB»WB,PB,NQ,D11 K 7 

COMMON CORD(100,2),NOP(200,4)»IMAT(200).ORT(25»5),NBC(25),NFIX(25) K 8 

l.Rl(200),R2(200).R3(200).R10(200),R20(200),R30(200).FORS(200),SM(2 K 9 

200, 15),SK(200, 15),ISP(200,15).SMPEM(200»15).ESTIF(12.12).EMASS(12» K 10 

312).NFIXK(25) K 11 

COMMON /COMP/ QBR(3,3,25)>.ABD(S.B),TH(25),ZK(25),MLAYER K 12 

DIMENSION 0(3,3), TK(25) K 13 
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DO 102 J=l»3 
DO 102 K=l,3 

ABD(J+3»K+3)=0.0 

ABD ( J+3 . K ) = ABD ( J+3 , K+3 ) 

ABD ( J , K+3 ) =ABD ( J+3 , K ) 

ABD(JfK)=ABD(J,K+3) 

DO 101 1=1,25 
DBR(J,K, I)=0. 

101 CONTINUE 

102 CONTINUE 

READ (5,108) MLAYER 
M=MLAYER 

READ (5,109) (L,TH(L),TK(L),I=1,M) 

TTK=0,0 
ZK(1)=TTK 
DO 103 1=1, M 
TTK=TTK+TK(I) 

2K(I+1)=TK(I)+ZK(I) 

103 CONTINUE 
MN=M+1 

DO 104 1=1, MN 

ZK(I)=ZK(I)-TTK/2. 

104 CONTINUE 
DEL=4.»ATAN(1. )/180. 

DEN= 1 . -ORT ( 1 , 2 ) *ORT ( 1 , 4 ) ««2/0RT C 1 , 1 ) 

Q(l,l)=ORT(l, D/DEN 
Q(2,2)=0RT(1,2)/DEN 
Q(2,1)-0RT(1,4)*Q(2,2) 

Q(1,2)=Q(2,1) 

Q(3,3)=0RT(1,3) 

Q(3,2)=0.0 

Q(3,1)=Q(3,2) 

Q(2,3)=Q(3,1) 

Q(1,3)=Q(2,3) 

DO 105 1=1, M 

ANGL=TH(I)»DEL 

C=COS(ANGL) 

W=SIN(ANGL) 

QBR(l,l,I)=Q(l,l)«C«»4+2,«(Qa,2)+2.»Q(3,3))»(C»W)«»2+a(2,2)*M* 
1 *4 

QBR(2,1,I)=(D(1, l)+Q(2,2)-4.»Q(3,3))*(C*W)**2+Q(l,2)*(W**4+C*»f4 
1 ) 

QBR(1,2,I)=0BR(2,1,I) 

QBR(2,E, I)=Q(1, 1)*W**4+2.*(Q(1,2)+2.»Q(3,3))»(C»W)»»2+Q(2,2)*C» 
1 »4 

QBR(3, 1, I)=(Q(l,l)-Q(l,2)-2.»Q(3,3))»W»C*^t3+(Q(l,2)-Q(2,2)+(2. ) 
1 *0(3,3))*(W)*(C**3) 

QBR(1,3, I)=DBR(3,1,I) 

QBR(3,2, 1) = (Q(1, l)-Q(l,8)-2.^fQ(3,3))*W«f*3*C+(Q(l,2)-Q(2,2)+2.'»Q 
1 (3,3))*W*C«»3 

QBR(E,3, I)=QBR(3,2,I) 

QBR(3,3,I)=(Q(1, l)+Q(2,2)-2.«a(l,2)-2,*-Q(3,3))*(W»C)**2+Q(3,3)» 
1 (W»»4+C*<t4) 

C 

C WRITE(B,500) I,TH(I),TK(I) 

C WRITE(B,510) 

C WRITE(B,5B0) ( (QBR( J, K, I ) , K=l, 3) , J=l,3) 

C 

105 CONTINUE 

DO 107 J=l,3 
DO 107 K=l,3 
DO lOG 1=1, M 

ABD ( J , K ) =ABD ( J , K ) +QBR ( J, K, I ) » ( ZK ( 1+ 1 ) -ZK ( I ) ) 
ABD(J,K+3)=ABD(J+3,K)+QBR(J,K,I)^f(ZK(I+l)^H^2-ZK(I)*-»2)/2. 

ABD ( J+3 , K) = ABD ( J , K+3 ) 

ABD ( J+3 , K+3 ) = ADD ( J+3 , K+3 ) +QBR ( J , K , I ) » ( ZK ( I + 1 ) »*3-ZK( I ) **3 ) /3 

1 

lOG CONTINUE 
107 CONTINUE 

WRITE (B,110) 
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NRITE 

(8,111) ((ABD(I,J),J=1,G),I=1,6) 

K 

84 




K 

85 


500 F0RMAT(2X.»LAYER=».I2,5X.*ftNGLE=»»F5.2.5X.*THICKNESS=«»F7.3) 

K 

86 


510 FDRMAT(2X,*QBAR-mTRIX») 

K 

87 


520 F0RMAT(5X»3E12.3/) 

K 

88 




K 

89 


RETURN 


K 

90 
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108 

FORMAT 

(15) 

K 

92 

109 

FORMAT 

(I5,F5.0,F10.0) 

K 

93 

no 

FORMAT 

(///,1X,10HABD MATRIX) 

K 

94 

111 

FORMAT 

(1X,BE11.3) 

K 

95 




K 
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APPENDIX B 

A COMPUTER PROGRAM FOR ESTIMATING THE CONTACT FORCE HISTORY BY USING THE 
EQUIVALENT MASS MODEL 

This program has been written for simply-supported beams only. 
Cantilever beams and simply-supported plates will be added in the 
near future. 

This program will be a subprogram in a large finite element 
program capable of analyzing impact responses of beams and plates. 

This subprogram will be used to provide an estimate of the contact 
time so that one may select a proper time increment for the finite 
difference used in the program. For this reason, the input cards 
for this subprogram were written to be identical to that for the 
program presented in Appendix A. 
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Listing of Program 

PROGRAM MAIN ( INPUT* OUTPUT, TAPE5=INPUT. TAPEB=OUTPUT) A 2 

COMMON /CQNTR/ TITLE(10),NP,NE.NB,NDF,NCN,NLD,NMAT,NSZF,LI,NT4,ND+ A 3 

1N,MATP,NPR0B A 4 

COMMON /TIME/ T.DT.DDT.TAU.KCON.KCNT A 5 

COMMON /DISP/ 01,02,03,010,020,030 A B 

COMMON /DIME/ TB,WB,PB,NQ,D11 A 7 

COMMON /SPHERE/ STF,R,CABU(10),QKONST(101 A 8 

COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25) A 9 

1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SMC2 A 10 

200,15),SK(200, 15),ISP(200, 15),SMPEM(200, 15),ESTIF(12, 12) , EMASSC 12, A 11 
312),NFIXK(25) A 12 

COMMON /COMP/ QBR(3,3,25),ABD(8,B),TH(25),2K(25),MLAYER A 13 

COMMON X1,X2,ND1,ND2 A 14 

REAL IB A 15 

A IB 

READ AND PRINT TITLE AND CONTROL A 17 

A 18 

101 READ (5,117) NPROB, (TITLE(I),I=1,10) A 19 

IF (NPROB. EO.O) GO TO 107 A 20 

NRITE (8,118) NPROB, (TITLE(I), 1=1, 10) A 21 

READ (5,108) NP,NE,NB,NTM,NMAT,NDIN,MATP,NDC,I1 A 22 

NDF=3 A 23 

NLD=NDIN*NTM A 24 

FLD=FLOAT(NLD)/10. A 25 

FDIN=FLOAT(NDIN)/10. A 28 

READ (5,114) TB,WB,R,N0,02,DT A 27 

WRITE (8,109) NP,NE,NB,NLD,NDF,NMAT A 28 

NLD=NLD+1 A 29 

A 30 

READ AND PRINT MATERIAL DATA A 31 

SPHERE DATA! L=NMAT (LAST MAT. CARD) A 32 

A 33 

READ (5,113) (L,(0RT(L,I),I=1,5),N=1,NMAT) A 34 

PB=0RT(NMAT,5) A 35 

WRITE (8,123) TB,WB,PB,R,NQ,Q2,DT A 38 

WRITE (8,122) A 37 

WRITE (8,118) (N, (0RT(N,I), 1=1,5), N=1,NMAT) A 38 

A 39 

READ INDENTATION CARD A 40 

A 41 

READ (5,112) STF A 42 

A 43 

READ NODAL POINT DATA A 44 

AND . A 45 

READ ELEMENT DATA A 48 

A 47 

DO 103 1=1, NDC A 48 

READ (5,115) ND1,ND2,X1,X2,IMT A 49 

EL=(X2-X1)/FL0AT(ND2-ND1) A 50 

C0RD(ND1,1)=X1 A 51 

C0RD(ND2, 1)=X2 A 52 

CORD(ND2,2)=0.0 A 53 

C0RD(ND1,2)=C0RD(ND2,2) A 54 

NDD=ND2-1 A 55 

DO 102 J=ND1,NDD A 58 

C0RD(J+1,1)=C0RD(J,1)+EL A 57 

CORD(J+1,2)=0.0 A 58 

NOP(J,l)=J A 59 

N0P(J,2)=J+1 A 80 

NOP(J,4)=0 A 81 

N0P(J,3)=N0P(J,4) A 82 

IMAT(J)=IMT A 83 

102 CONTINUE A 84 

103 CONTINUE A 85 

A 86 

READ BOUNDARY DATA A 67 

A 68 

READ (5,111) (NBC(I),NFIX(I),I=1,NB)' A 69 

IF (MATP.EQ.l) CALL CMPD A 70 

C A 71 
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ISOTROPIC MATP=0.0 A 

COMPOSITE MATP=1.0 A 

A 

IF (Il.NE.O) GO TO 104 A 

A 

PRINT INPUT DATA A 

A 

WRITE (6,119) A 

WRITE (6,110) (N, (C0RD(N,M),M=1,2),N=1,NP) A 

A 

WRITE(6,103) A 

WRITE(6,3)(N, (N0P(N,M),M=1,4),IMAT(N),N=1,NE) A 

A 

WRITE (6,120) A 

WRITE (6,111) (NBC(I),NFIX(I),I=1.NB) A 

WRITE (6,121) FDIN,FLD A 

104 CONTINUE A 

DO 105 IJ=1,200 A 

R10(IJ)=0. A 

R20(IJ)=0. A 

R30(IJ)=0. A 

105 FORS(IJ)=0. A 

DO 106 IJ=1,25 A 

lOB NFIXK(IJ)=NFIX(IJ) A 

CALL TMX A 

A 

3 FORMAT (615) A 

A 

103 FORMATdOHO ELEMENTS/9X, 1HI,4X, 1HJ,4X, 1HK,8X,3HMAT) A 

A 

GO TO 101 A 

107 STOP A 

C A 

108 FORMAT 015) A 

109 FORMAT (13HONODAL POINTS, 9X, I5/1X,8HELEMENTS, 13X, I5/1X, 19IHB0UNDARY A 

1 CONDITIONS, 2X,I5/1X,18H0UTPUT LIMIT, lOX, I5/1X, 18HDEGREES OF FREED A 
20M,3X,I5/1X,9HMATERIALS, 12X,I5) A 

110 FORMAT (I10,2F10.3) A 

111 FORMAT (215) A 

112 FORMAT (E10.3) A 

113 FORMAT (I5,5F10.0) A 

114 FORMAT (3F10.0/I5,2F10.0) A 

115 FORMAT (215, 2F10. 0,15) A 

116 FORMAT (I5,7X,3(F10.1,4X),F5.3,7X,F8.B//) A 

117 FORMAT (12, 10A7) A 

118 FORMAT (1H1,I2,10A7) A 

119 FORMAT (14H0 NODAL POINTS/ITX, IHX, lOX, IHY) A 

120 FORMAT (21H0 BOUNDARY CONDITIONS) A 

121 FORMAT (ISHOPRINTlNG SCHEME/5X,22H1. REPORT OUTPUT EUERY,F6.2,2X,4 A 

1HMSEC/5X.22H2. TERMINATE OUTPUT AT, F6.2, 2X, 4HMSEC) A 

122 FORMAT (1H0.20H MATERIAL PROPERTIES/IX, 8HMAT. N0.,7X,2HE1, 12X,2HE2 A 

1, 11X.3HG12, i0X,3HU12, 10X,3HRHO/) A 

123 FORMAT (19H0BEAM THICKNESS(IN),7X,FB.3/1X, 14HBEAM WIDTH! IN) , 1 IX, FB A 

1.3/1X,22HSPHEREDENSITY(SL/IN3),4X,F8.B/1X,17HSPHERERADIUS(IN),8X A 
2,F6.3//1X, IIHIMPACT NODE, 14X, I2/1X, 23H1MPACT UELOCITY! IN/SEC) , 2X, F A 
36.1/lX, 39H INTEGRA! I ON TIME INCREMENT! X E-OB SEC),E10.3) A 

C A 

END A 

SUBROUTINE CMPD B 

COMMON /CONTR/ TITLE(10),NP,NE,NB,NDF,NCN,NLD,NMAT,NSZF,LI,NT4,NDI B 
1N,MATP,NPR0B B 

COMMON /TIME/ T, DT, DDT, TAU, KCON, KCNT B 

COMMON /DISP/ 01, 02, 03, 010,020,030 B 

COMMON /DIME/ TB,WB,PB,NQ,D11 B 

COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25) B 
1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 B 
200, 15),SK(200, 15),ISP(200, 15),SMPEM(200, 15), ESTIF( 12, 12) , EMASS( 12, B 
312),NFIXK(25) B 

COMMON /COMP/ QBR(3,3,25),ABD(6,6),TH(25),ZK(25),MLAYER B 

COMMON X1,X2,ND1,ND2 B 
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REAL IB 

DIMENSIOM Q(3i3). TK(25) 

DO 102 J=l,3 
DO 102 K=l,3 

ABD(J+3.K+3)=0.0 

ABD ( J+3 , K ) = ABD ( J+3. K+3 > 

ABD(J.K+3)=ABD(J+3,K) 

ABD(J.K)=ABD(J,K+3) 

DO 101 1=1,25 
QBR(J,K,I)=0. 

101 CONTIMUE 

102 CONTIHUE 

READ (5,108) MLAYER 
t1=MLAYER 

READ (5,109) (L,TH(L),TK(L),I=1,M) 

TTK=0.0 
ZK(1)=TTK 
DO 103 1=1, M 
TTK=TTK+TK(I) 

ZK(I+1)=TK(I)+ZK(I) 

103 COMTINUE 
MM=M+1 

DO 104 1=1, HM 

ZK(I)=ZK(l)-TTK/2. 

104 CONTINUE 
DEL=4.*ATAN(1. )/180. 

DEN= 1 . -ORT ( 1 , 2 ) *ORT ( 1 , 4 ) *»2/0RT ( 1 , 1 ) 

0(1, l)=ORT(l, D/DEN 
Q(2,2)=0RT(1,2)/DEN 
0(2, 1)=0RT(1,4)*Q(2,2) 

0(1,2)=Q(2,1) 

0(3,3)=0RT(1,3) 

0(3,2)=0.0 
0(3, 1)=0(3,2) 

0(2,3)=Q(3, 1) 

0(1,3)=0(2,3) 

DO 105 1=1, M 

ANGL=TH(I)*DEL 

C=COS(ANGL) 

U=SIN(ANGL) 

QBR( 1 , 1 , I )=Q( 1, 1 )*C*»4+2.»(Q( 1,2)+2.*Q(3,3) )»(C*W)*»2+Q(2,2)»W» 
1 *4 

QBR(2,l,I)=(Q(l,l)+0(2,2)-4.»Q(3,3))*(C*U)»»2+Q(l,2)»(W**4+C»*4 
1 ) 

QBR(1,2,I)=0BR(2,1,I) 

QBR(2,2, 1)=Q(1,1)^«■^I**4+2.*(Q(1,2)+2.*0(3,3))»(C*^I)**2+Q(2,^)*C» 
1 *4 

OBR(3, 1, I)=(Q(l,l)-Q(l,2)-2,^»Q(3,3))*W»C**3+(Q(l,2)-Q(2,2)+(2.) 
1 *a(3,3))*(W).*(C«*3) 

QBR(1,3, I)=0BR(3, 1,1) 

0BR(3,2, I)=(Q(1, 1)-0(1,2)-2.*Q(3,3))*W*«-3*C+(0(1,2)-Q(2,2)+2.*Q 
1 (3,3))*W*C«»3 

QBR(2,3,I)=QBR(3,2,I) 

0BR(3,3, I)=(0(1, l)+0(2,2)-2.»Q(l,2)-2.*Q(3,3))*(N»C)*f*2+D(3,3)* 
1 (W*M+C**4) 


MRITE(6,500) I,TH(I),TK(I) 

MRITE(G,510) 

WRITE(G,520) ( (OBR( J, K, I ), K=l,3) , J=l, 3) 


105 CONTINUE 

DO 107 J=l,3 
DO 107 K=l,3 
DO lOG 1=1, N 

ABD(J,K)=ABD(J,K)+OBR(J,K,I)«(ZK(I+l)-ZK(I)) 

ABD(J,K+3)=ABD(J+3,K)+QBR(J,K,I)»(ZK(I+l)»*2-ZK(I)»*2)/2.. 

ABD(J+3,K)=ABD(J,K+3) 

ABD(J+3,K+3)=ABD(J+3,K+3)+QBR(J,K, I)»(ZK(I+l)»*3-ZK(I)**3)/3 
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CONTINUE 
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102 


107 COMTIISUE B 

WRITE ((ABDa,J),J=l.G)*I=l,B) B 

WRITE B 

B 

500 F0RMAT(2Xf ^»LAYER=». 12, 5X, »AMGLE=*, F5. 2, 5X, *THICKMESS=«, F7. 3) B 
510 F0RMAT(2X,*QBAR-MATRIX») B 

520 FDRMAT(5X,3E12.3/) B 

B 

RETURH B 

B 

108 FORMAT (15) B 

109 FORMAT a5,F5.Q,F10.0) B 

110 FORMAT (///, IX, lOHABD MATRIX) B 

111 FORMAT (1X,GE11.3) B 

B 

EMD B 

SUBROUTINE TMX C 

COMMON /CONTR/ TITLE(10),NP,NE,NB,NDF,NCN,NLD,NMAT,NSZF,LI,NT4,NDI C 
1N,MATP,NPR0B C 

COMMON /TIME/ T,DT,DDT,TAU,KCON,KCNT C 

COMMON /DISP/ Q1,Q2, 03, 010,020,030 C 

COMMON /BIMB/ TB, WB, PB, NO, Dll C 

COMMON /SPHERE/ STF,R,CABU(10),QKONST(10) C 

COMMON /PLASTIC/ DISPEM,NDISPEM,FORSPM,DISPM C 

COMMON CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25) C 

1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 C 
200, 15),SK(200,15),ISP(200,15),SMPEM(200, 15) , ESTIFC 12, 12),EMASS(12, C 
312),NFIXK(25) C 

COMMON /COMP/ QBR(3,3,25),ABD(B,B),TH(25),ZK(25),MLAYER C 

COMMON X1,X2,ND1,ND2 C 

REAL IB C 

C 

EQUIUALENT MODEL FOR ESTIMATING TIME INTERUAL C 

C 

IF (STF.NE.0.0) GO TO 101 C 

STFI=(4./3.)*SQRT(R)/((1.-ORT(NMAT,4)»*2)/ORT(NMAT,1)+(1.-ORT(1,4) C 
1*»2)/0RT(1, D) C 

gTFA=(4./3.)»SQRT(R)/((l.-0RT(NMAT,4)»»2)/DRT(NMAT,l)+l./0RT(l,2)) C 
STF=STFI C 

IF (MATP.EQ.l) STF=STFA C 

101 PAI=4.*ATAN(1. ) C 

B ALLM= ( 4 . /3 . ) »P A I * ( R«*3 ) *PB C 

BL=X2-X1 C 

AB=WB»TB C 

IB=WB*TB*«3/12. C 

WATP=FLOAT(MATP) C 

Dll=ORT(l, 1)*IB C 

IF (MATP.EQ.l) D11=ABD(4,4) C 

WRITE (G,105) BL,STF,WATP,D11,ABD(4,4) C 

WN 1= ( (Dll * ( PAI /BL ) **4 ) / ( AB»ORT (1,5)))»»0.5 C 

TN1=2.*PAI/WN1 C 

TD=0.0 C 

DO 103 1=1,1000 C 

N=0 C 

TD=TD+1.0E-7 C 

F1=0.0 C 

G1=0.0 C 

DO 102 J=l,100 C 

N=N+1 C 

C=FL0AT(NQ-1)/FL0AT(ND2-ND1) C 

WX=SIN(N*PAI*C) C 

PP=1./(4.»N'»»4*TD**2-TN1»»2) C 

QQ=TN1/(2.<»N»«-2*TD) C 

SS=WNl*fN«»»2»TD/2. C 

SUMF=(PP*N^J*2*(1-QQ*SIN(SS))»WX)»*2 C 

SUMG=(PP«C0S(SS)*WX)**2 C 

F1=F1+SUMF C 

G1=G1+SUMG C 

102 CONTINUE C 

SUM 1=2 . ^^F 1 ■» IG . 3«TD*-*2/ ( D 1 1*PAI »»2 ) C 
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103 


SUM2=2.*Gl<flS.*0RT(l,5)*ftB*BL**7ADll*»2*PPlI**4) C 5G 

EMT=l./(SUm+SUM2) C 57 

STFE= 1 . /EMT+ 1 . /BALLM C 58 

STFT=STF*STFE C 59 

APHAmX=(1.25*Q2**2/STFT)**0.4 C BO 

T0TALT=2.94*APHAMAX/Q2 C B1 

FMAX=STF*APHAMAX**1.5 C B2 

EPS=TOTALT-TD C B3 

IF (ABS(EPS).LE.l.E-?) GO TO 104 C B4 

103 CONTINUE C B5 

104 WRITE (S,10B) TOTALT,TD,FNAX» APHAMAX*ENT C BB 

RETURN C B7 

C C 68 

105 FORMAT (//,5X, 5HBEAM=, E10.3. 5X. 4HSTF=,E10.3»5X. 5HMATP=.E10.3 C B9 

1»5X» 4HD11=,E10.3.5X, 9HABD(4,4)=.E10.3) C 70 

lOB FORMAT C//,5X, 7HTDTALT=,E10.3,5X, 3HTD=,E10.3.5Xt 5HFMAX=,E10. C 71 

13.5X, 8HAPHAMAX=,E10.3,5X. 4HEMT=,E10.3) C 72 

C C 73 

END C 74 
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